package 'xlsx' successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\cvasquezv\AppData\Local\Temp\Rtmpm493xe\downloaded_packages
Experimentos repetidos en el espacio
Análisis de DCA de dos vías repetido en el espacio
modelo.dca1 <-lm(Diametro_tallo_90dds ~ Zona + BLOQUE %in% Zona + VARIEDAD + VARIEDAD/Zona, data = data_D)modelo.dca2 <-lm(Diametro_tallo_90dds ~ Zona + BLOQUE %in% Zona + VARIEDAD, data = data_D)
Nota
La expresión “BLOQUE/ZONA” es para considerar la interacción como efecto fijo.
La expresión “BLOQUE:ZONA” es para considerar la interacción como efecto aleatorio.
La expresión “BLOQUE %in% ZONA” permite que se genere una estructura del análisis de varianza donde se anide el efecto de la variedad dentro de cada bloque
Durbin-Watson test
data: modelo.dca
DW = 1.8537, p-value = 0.9119
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
shapiro.test(rstudent(modelo.dca))
Shapiro-Wilk normality test
data: rstudent(modelo.dca)
W = 0.98058, p-value = 0.9859
ad.test(rstudent(modelo.dca))
Anderson-Darling normality test
data: rstudent(modelo.dca)
A = 0.19482, p-value = 0.861
lillie.test(rstudent(modelo.dca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dca)
D = 0.12277, p-value = 0.8862
Exact one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dca)
D = 0.11883, p-value = 0.9881
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dca))
Cramer-von Mises normality test
data: rstudent(modelo.dca)
W = 0.029639, p-value = 0.8354
pearson.test(rstudent(modelo.dca))
Pearson chi-square normality test
data: rstudent(modelo.dca)
P = 2, p-value = 0.5724
sf.test(rstudent(modelo.dca))
Shapiro-Francia normality test
data: rstudent(modelo.dca)
W = 0.97048, p-value = 0.8503
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 1.889962, Df = 1, p = 0.16921
bptest(modelo.dca)
studentized Breusch-Pagan test
data: modelo.dca
BP = 12, df = 6, p-value = 0.06197
bptest(modelo.dca, studentize = F)
Breusch-Pagan test
data: modelo.dca
BP = 5.9393, df = 6, p-value = 0.43
olsrr::ols_test_breusch_pagan(modelo.dca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
------------------------------------------------
Response : Diametro_tallo_90dds
Variables: fitted values of Diametro_tallo_90dds
Test Summary
----------------------------
DF = 1
Chi2 = 1.889962
Prob > Chi2 = 0.1692062
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.
Estadísticas globales
modelo.dca %>%gvlma()
Call:
lm(formula = Diametro_tallo_90dds ~ Zona + BLOQUE %in% Zona +
VARIEDAD, data = data_D)
Coefficients:
(Intercept) ZonaCastilla VARIEDAD2
0.72083 -0.07500 -0.04167
ZonaCamana:BLOQUEII ZonaCastilla:BLOQUEII ZonaCamana:BLOQUEIII
0.05000 0.17500 0.07500
ZonaCastilla:BLOQUEIII
0.05000
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = .)
Value p-value Decision
Global Stat 5.213e+00 0.2661 Assumptions acceptable.
Skewness 1.568e-32 1.0000 Assumptions acceptable.
Kurtosis 5.102e-01 0.4751 Assumptions acceptable.
Link Function 2.589e+00 0.1076 Assumptions acceptable.
Heteroscedasticity 2.113e+00 0.1460 Assumptions acceptable.
\(Y_{ijk}\) = Valor observado de la variable respuesta.
\(\hat{Y}_{ijk}\) = Valor ajustado de la variable respuesta.
\(\mu\) = Promedio observado de la variable respuesta.
\(\tau_{i}\) = Efecto del i-ésimo nivel del factor Zona.
\(\text{Error}(\tau\text{rep})_{i(k)}\) = Efecto del i-ésimo nivel de Zona en el k-ésimo nivel de repetición.
\(\beta_{j}\) = Efecto del j-ésimo nivel del factor tratamiento (Variedad).
\(\epsilon_{ijk}\) = Residuo observado del modelo.
Pruebas de hipótesis
Para el factor A (Zona):
\(H_0: \tau_{A1} = \tau_{A2} = 0\)
\(H_1: \text{En al menos un nivel del factor A el } \tau \text{ es diferente a los demás.}\)
\(H_1: \tau_i \neq 0\text{; en al menos un nivel del factor A.}\)
Para el factor B (Variedad):
\(H_0: \beta_{B1} = \beta_{B2} = 0\)
\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)
\(H_1: \beta_j \neq 0\text{; en al menos un nivel del factor B.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
anova(modelo.dca, test ="F")
Analysis of Variance Table
Response: Diametro_tallo_90dds
Df Sum Sq Mean Sq F value Pr(>F)
Zona 1 0.005208 0.0052083 3.0488 0.14123
VARIEDAD 1 0.005208 0.0052083 3.0488 0.14123
Zona:BLOQUE 4 0.038333 0.0095833 5.6098 0.04315 *
Residuals 5 0.008542 0.0017083
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
GAD::estimates(modelo.dca)
$tm
Zona BLOQUE VARIEDAD n
Zona 0 3 2 1.714286
VARIEDAD 2 3 0 1.714286
Zona:BLOQUE 1 1 2 1.714286
Res 1 1 1 1.000000
$mse
Mean square estimates
Zona "Res + Zona:BLOQUE + Zona"
VARIEDAD "Res + VARIEDAD"
Zona:BLOQUE "Res + Zona:BLOQUE"
Residual "Res"
$f.versus
F-ratio versus
Zona "Zona:BLOQUE"
VARIEDAD "Residual"
Zona:BLOQUE "Residual"
GAD::gad(modelo.dca)
Analysis of Variance Table
Response: Diametro_tallo_90dds
Df Sum Sq Mean Sq F value Pr(>F)
Zona 1 0.005208 0.0052083 0.5435 0.50190
VARIEDAD 1 0.005208 0.0052083 3.0488 0.14123
Zona:BLOQUE 4 0.038333 0.0095833 5.6098 0.04315 *
Residual 5 0.008542 0.0017083
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Importante
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(Zona/rep)” en el caso de efectos aleatorios o “Error(Zona:rep)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos.
aov(Diametro_tallo_90dds ~ Zona + VARIEDAD +Error(BLOQUE %in% Zona), data = data_D) -> aov.dcasummary(aov.dca)
Error: BLOQUE:Zona
Df Sum Sq Mean Sq F value Pr(>F)
Zona 1 0.00521 0.005208 0.543 0.502
Residuals 4 0.03833 0.009583
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
VARIEDAD 1 0.005208 0.005208 3.049 0.141
Residuals 5 0.008542 0.001708
broom::tidy(aov.dca)
# A tibble: 4 × 7
stratum term df sumsq meansq statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 BLOQUE:Zona Zona 1 0.00521 0.00521 0.543 0.502
2 BLOQUE:Zona Residuals 4 0.0383 0.00958 NA NA
3 Within VARIEDAD 1 0.00521 0.00521 3.05 0.141
4 Within Residuals 5 0.00854 0.00171 NA NA
broom::tidy(gad(modelo.dca))
# A tibble: 4 × 6
term df sumsq meansq statistic p.value
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 Zona 1 0.00521 0.00521 0.543 0.502
2 VARIEDAD 1 0.00521 0.00521 3.05 0.141
3 Zona:BLOQUE 4 0.0383 0.00958 5.61 0.0432
4 Residual 5 0.00854 0.00171 NA NA
Valor de la tabla de F para el factor Zona con una significancia de 0.05.
qf(0.95, 1, 4)
[1] 7.708647
Valor de la tabla de F para el factor Variedad con una significancia de 0.05.
qf(0.95, 1, 5)
[1] 6.607891
Conclusión.
Con respecto al Factor Zona: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor Zona tienen un efecto sobre el diámetro de tallo estadísticamente similar a 0.
Con respecto al Factor Variedad: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor Variedad tienen un efecto sobre el diámetro de tallo estadísticamente similar a 0.
agricolae::cv.model(modelo.dca)
[1] 5.733918
Comparaciones de medias para los efectos principales del Factor Zona
data %>%with(LSD.test( Diametro_tallo_90dds, # Cambiar según nombre de variable respuesta Zona, # Cambiar según nombre de variable independienteDFerror =get_df_ea(aov.dca), MSerror =get_mse_ea(aov.dca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: Diametro_tallo_90dds ~ Zona
LSD t Test for Diametro_tallo_90dds
Mean Square Error: 0.009583333
Zona, means and individual ( 95 %) CI
Diametro_tallo_90dds std r LCL UCL Min Max
Camana 0.7270833 0.09205288 24 0.6716027 0.782564 0.5 0.90
Castilla 0.7395833 0.06753287 24 0.6841027 0.795064 0.6 0.85
Alpha: 0.05 ; DF Error: 4
Critical Value of t: 2.776445
least Significant Difference: 0.07846153
Treatments with the same letter are not significantly different.
Diametro_tallo_90dds groups
Castilla 0.7395833 a
Camana 0.7270833 a
Comparaciones de medias para los efectos principales del Factor Variedad
data %>%with(LSD.test( Diametro_tallo_90dds, # Cambiar según nombre de variable respuesta VARIEDAD, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dca), MSerror = dvmisc::get_mse(modelo.dca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: Diametro_tallo_90dds ~ VARIEDAD
LSD t Test for Diametro_tallo_90dds
Mean Square Error: 0.001708333
VARIEDAD, means and individual ( 95 %) CI
Diametro_tallo_90dds std r LCL UCL Min Max
1 0.7583333 0.07322786 24 0.7366457 0.780021 0.6 0.90
2 0.7083333 0.08030738 24 0.6866457 0.730021 0.5 0.85
Alpha: 0.05 ; DF Error: 5
Critical Value of t: 2.570582
least Significant Difference: 0.03067094
Treatments with the same letter are not significantly different.
Diametro_tallo_90dds groups
1 0.7583333 a
2 0.7083333 b
Comparaciones de medias para las interacciones
Para los niveles del factor A dentro del nivel B1:
A1 vs A2:
\(H_0: \mu_{A1} - \mu_{A2} = 0\)
\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)
Para los niveles del factor B dentro del nivel A1:
B1 vs B2:
\(H_0: \mu_{B1} - \mu_{B2} = 0\)
\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)
NOTA: Repetir este proceso para cada nivel de A y cada nivel de B.
Análisis de varianza para interacción de dos factores con el paquete phia
Comparación de los niveles de B dentro de cada nivel de A
multcomp.test_2factors <-function(object, respuesta, factor_name1, factor_name2, test, aov){# Función auxiliar para aplicar la prueba de comparaciones múltiples a un data frame multcomp_df <-function(df, respuesta, factor_name, test, aov){if (test =="LSD") { comp <-LSD.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[5]] } elseif (test =="HSD") { comp <-HSD.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[5]] } elseif (test =="duncan") { comp <-duncan.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[6]] } elseif (test =="SNK") { comp <-SNK.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[6]] } else {stop("Test no válido. Los tests disponibles son LSD, HSD, duncan, y SNK.") }return(comp %>%rename("y"="df[[respuesta]]") ) }# Aplicar la prueba de comparaciones múltiples a todos los data frames dentro de filters1 comp_filters1 <-lapply(object[[1]], function(df){multcomp_df(df, respuesta, factor_name2, test, aov) %>%arrange(row.names(.)) %>%rownames_to_column(var ="x") %>%relocate(x) })# Aplicar la prueba de comparaciones múltiples a todos los data frames dentro de filters2 comp_filters2 <-lapply(object[[2]], function(df){multcomp_df(df, respuesta, factor_name1, test, aov) %>% dplyr::arrange(row.names(.)) %>%rownames_to_column(var ="x") %>%relocate(x) %>% dplyr::mutate(groups =toupper(groups)) })row.names(comp_filters1) <-NULLrow.names(comp_filters2) <-NULL# Retornar una lista con las comparaciones múltiples para cada data frame result <-list() result[[as.name(substitute(factor_name1))]] <- comp_filters1 result[[as.name(substitute(factor_name2))]] <- comp_filters2return(#list(#"Comparación de los niveles del factor B dentro de cada nivel del factor A", result#[[1]],# "Comparación de los niveles del factor A dentro de cada nivel del factor B",# result[[2]]# ) )}
$Zona
$Zona$Camana
x y groups
1 1 0.7708333 a
2 2 0.6833333 b
$Zona$Castilla
x y groups
1 1 0.7458333 a
2 2 0.7333333 a
$VARIEDAD
$VARIEDAD$`1`
x y groups
1 Camana 0.7708333 A
2 Castilla 0.7458333 A
$VARIEDAD$`2`
x y groups
1 Camana 0.6833333 B
2 Castilla 0.7333333 A
# Convertir la lista a un data framedf <- result.comp %>% reshape2::melt() %>%ungroup() %>% dplyr::select(-variable) %>%relocate("Factor"="L1","Nivel"="L2", x,"y"="value", groups)df
Factor Nivel x y groups
1 Zona Camana 1 0.7708333 a
2 Zona Camana 2 0.6833333 b
3 Zona Castilla 1 0.7458333 a
4 Zona Castilla 2 0.7333333 a
5 VARIEDAD 1 Camana 0.7708333 A
6 VARIEDAD 1 Castilla 0.7458333 A
7 VARIEDAD 2 Camana 0.6833333 B
8 VARIEDAD 2 Castilla 0.7333333 A
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# # Convertir level1 y level2 en nombres simbólicos# level1 <- as.name(level1)# level2 <- as.name(level2)# # filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1, level2)) %>% dplyr::mutate(groups =paste0(groups.y,groups.x)) %>% dplyr::select(!!level1,!!level2, y, groups) #%>%# rename(!!level1 := x,# !!level2 := Nivel)# devolver una lista con los dos subconjuntosreturn(df)}
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1,level2)) %>% dplyr::mutate(y =paste0(round(y,2), " ", groups.y, groups.x)) %>% dplyr::select(!!level1,!!level2, y)# devolver una lista con los dos subconjuntosreturn(df)}
Análisis de DBCA de dos vías repetido en el espacio
Definición del modelo
modelo.dbca1 <-lm(Diametro_tallo_90dds ~ BLOQUE + Zona + BLOQUE %in% Zona + VARIEDAD + BLOQUE + VARIEDAD/Zona, data = data_D)modelo.dbca2 <-lm(Diametro_tallo_90dds ~ Zona + BLOQUE %in% Zona + VARIEDAD, data = data_D)
Nota
La expresión “BLOQUE/ZONA” es para considerar la interacción como efecto fijo.
La expresión “BLOQUE:ZONA” es para considerar la interacción como efecto aleatorio.
La expresión “BLOQUE %in% ZONA” permite que se genere una estructura del análisis de varianza donde se anide el efecto de la variedad dentro de cada bloque
Durbin-Watson test
data: modelo.dbca
DW = 1.8537, p-value = 0.9119
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
shapiro.test(rstudent(modelo.dbca))
Shapiro-Wilk normality test
data: rstudent(modelo.dbca)
W = 0.98058, p-value = 0.9859
ad.test(rstudent(modelo.dbca))
Anderson-Darling normality test
data: rstudent(modelo.dbca)
A = 0.19482, p-value = 0.861
lillie.test(rstudent(modelo.dbca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dbca)
D = 0.12277, p-value = 0.8862
Exact one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dbca)
D = 0.11883, p-value = 0.9881
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dbca))
Cramer-von Mises normality test
data: rstudent(modelo.dbca)
W = 0.029639, p-value = 0.8354
pearson.test(rstudent(modelo.dbca))
Pearson chi-square normality test
data: rstudent(modelo.dbca)
P = 2, p-value = 0.5724
sf.test(rstudent(modelo.dbca))
Shapiro-Francia normality test
data: rstudent(modelo.dbca)
W = 0.97048, p-value = 0.8503
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dbca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 1.889962, Df = 1, p = 0.16921
bptest(modelo.dbca)
studentized Breusch-Pagan test
data: modelo.dbca
BP = 12, df = 6, p-value = 0.06197
bptest(modelo.dbca, studentize = F)
Breusch-Pagan test
data: modelo.dbca
BP = 5.9393, df = 6, p-value = 0.43
olsrr::ols_test_breusch_pagan(modelo.dbca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
------------------------------------------------
Response : Diametro_tallo_90dds
Variables: fitted values of Diametro_tallo_90dds
Test Summary
----------------------------
DF = 1
Chi2 = 1.889962
Prob > Chi2 = 0.1692062
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.
Estadísticas globales
modelo.dbca %>%gvlma()
Call:
lm(formula = Diametro_tallo_90dds ~ Zona + BLOQUE %in% Zona +
VARIEDAD, data = data_D)
Coefficients:
(Intercept) ZonaCastilla VARIEDAD2
0.72083 -0.07500 -0.04167
ZonaCamana:BLOQUEII ZonaCastilla:BLOQUEII ZonaCamana:BLOQUEIII
0.05000 0.17500 0.07500
ZonaCastilla:BLOQUEIII
0.05000
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = .)
Value p-value Decision
Global Stat 5.213e+00 0.2661 Assumptions acceptable.
Skewness 1.568e-32 1.0000 Assumptions acceptable.
Kurtosis 5.102e-01 0.4751 Assumptions acceptable.
Link Function 2.589e+00 0.1076 Assumptions acceptable.
Heteroscedasticity 2.113e+00 0.1460 Assumptions acceptable.
\(Y_{ijk}\) = Valor observado de la variable respuesta.
\(\hat{Y}_{ijk}\) = Valor ajustado de la variable respuesta.
\(\mu\) = Promedio observado de la variable respuesta.
\(\tau_{i}\) = Efecto del i-ésimo nivel del factor Zona.
\(\text{Error}(\tau\text{Bloque})_{i(k)}\) = Efecto del i-ésimo nivel de Zona en el k-ésimo nivel de repetición.
\(\beta_{j}\) = Efecto del j-ésimo nivel del factor tratamiento (Variedad).
\(\gamma_{k}\) = Efecto del k-ésimo nivel del factor Bloque.
\(\epsilon_{ijk}\) = Residuo observado del modelo.
Pruebas de hipótesis
Para el factor A (Zona):
\(H_0: \tau_{A1} = \tau_{A2} = 0\)
\(H_1: \text{En al menos un nivel del factor A el } \tau \text{ es diferente a los demás.}\)
\(H_1: \tau_i \neq 0\text{; en al menos un nivel del factor A.}\)
Para el factor B (Variedad):
\(H_0: \beta_{B1} = \beta_{B2} = 0\)
\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)
\(H_1: \beta_j \neq 0\text{; en al menos un nivel del factor B.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
anova(modelo.dbca, test ="F")
Analysis of Variance Table
Response: Diametro_tallo_90dds
Df Sum Sq Mean Sq F value Pr(>F)
Zona 1 0.005208 0.0052083 3.0488 0.14123
VARIEDAD 1 0.005208 0.0052083 3.0488 0.14123
Zona:BLOQUE 4 0.038333 0.0095833 5.6098 0.04315 *
Residuals 5 0.008542 0.0017083
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
GAD::estimates(modelo.dbca)
$tm
Zona BLOQUE VARIEDAD n
Zona 0 3 2 1.714286
VARIEDAD 2 3 0 1.714286
Zona:BLOQUE 1 1 2 1.714286
Res 1 1 1 1.000000
$mse
Mean square estimates
Zona "Res + Zona:BLOQUE + Zona"
VARIEDAD "Res + VARIEDAD"
Zona:BLOQUE "Res + Zona:BLOQUE"
Residual "Res"
$f.versus
F-ratio versus
Zona "Zona:BLOQUE"
VARIEDAD "Residual"
Zona:BLOQUE "Residual"
GAD::gad(modelo.dbca)
Analysis of Variance Table
Response: Diametro_tallo_90dds
Df Sum Sq Mean Sq F value Pr(>F)
Zona 1 0.005208 0.0052083 0.5435 0.50190
VARIEDAD 1 0.005208 0.0052083 3.0488 0.14123
Zona:BLOQUE 4 0.038333 0.0095833 5.6098 0.04315 *
Residual 5 0.008542 0.0017083
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Importante
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(Zona/rep)” en el caso de efectos aleatorios o “Error(Zona:rep)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos.
aov(Diametro_tallo_90dds ~ BLOQUE + Zona + VARIEDAD +Error(BLOQUE %in% Zona), data = data_D) -> aov.dbcasummary(aov.dbca)
Error: BLOQUE:Zona
Df Sum Sq Mean Sq F value Pr(>F)
BLOQUE 2 0.025417 0.012708 1.968 0.337
Zona 1 0.005208 0.005208 0.806 0.464
Residuals 2 0.012917 0.006458
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
VARIEDAD 1 0.005208 0.005208 3.049 0.141
Residuals 5 0.008542 0.001708
broom::tidy(aov.dbca)
# A tibble: 5 × 7
stratum term df sumsq meansq statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 BLOQUE:Zona BLOQUE 2 0.0254 0.0127 1.97 0.337
2 BLOQUE:Zona Zona 1 0.00521 0.00521 0.806 0.464
3 BLOQUE:Zona Residuals 2 0.0129 0.00646 NA NA
4 Within VARIEDAD 1 0.00521 0.00521 3.05 0.141
5 Within Residuals 5 0.00854 0.00171 NA NA
broom::tidy(gad(modelo.dbca))
# A tibble: 4 × 6
term df sumsq meansq statistic p.value
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 Zona 1 0.00521 0.00521 0.543 0.502
2 VARIEDAD 1 0.00521 0.00521 3.05 0.141
3 Zona:BLOQUE 4 0.0383 0.00958 5.61 0.0432
4 Residual 5 0.00854 0.00171 NA NA
Valor de la tabla de F para el factor Zona con una significancia de 0.05.
qf(0.95, 1, 2)
[1] 18.51282
Valor de la tabla de F para el factor Variedad con una significancia de 0.05.
qf(0.95, 1, 5)
[1] 6.607891
Conclusión.
Con respecto al Factor Zona: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor Zona tienen un efecto sobre el diámetro de tallo estadísticamente similar a 0.
Con respecto al Factor Variedad: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor Variedad tienen un efecto sobre el diámetro de tallo estadísticamente similar a 0.
agricolae::cv.model(modelo.dbca)
[1] 5.733918
Comparaciones de medias para los efectos principales del Factor Zona
data %>%with(LSD.test( Diametro_tallo_90dds, # Cambiar según nombre de variable respuesta Zona, # Cambiar según nombre de variable independienteDFerror =get_df_ea(aov.dbca), MSerror =get_mse_ea(aov.dbca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: Diametro_tallo_90dds ~ Zona
LSD t Test for Diametro_tallo_90dds
Mean Square Error: 0.006458333
Zona, means and individual ( 95 %) CI
Diametro_tallo_90dds std r LCL UCL Min Max
Camana 0.7270833 0.09205288 24 0.6565018 0.7976648 0.5 0.90
Castilla 0.7395833 0.06753287 24 0.6690018 0.8101648 0.6 0.85
Alpha: 0.05 ; DF Error: 2
Critical Value of t: 4.302653
least Significant Difference: 0.09981732
Treatments with the same letter are not significantly different.
Diametro_tallo_90dds groups
Castilla 0.7395833 a
Camana 0.7270833 a
Comparaciones de medias para los efectos principales del Factor Variedad
data %>%with(LSD.test( Diametro_tallo_90dds, # Cambiar según nombre de variable respuesta VARIEDAD, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dbca), MSerror = dvmisc::get_mse(modelo.dbca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: Diametro_tallo_90dds ~ VARIEDAD
LSD t Test for Diametro_tallo_90dds
Mean Square Error: 0.001708333
VARIEDAD, means and individual ( 95 %) CI
Diametro_tallo_90dds std r LCL UCL Min Max
1 0.7583333 0.07322786 24 0.7366457 0.780021 0.6 0.90
2 0.7083333 0.08030738 24 0.6866457 0.730021 0.5 0.85
Alpha: 0.05 ; DF Error: 5
Critical Value of t: 2.570582
least Significant Difference: 0.03067094
Treatments with the same letter are not significantly different.
Diametro_tallo_90dds groups
1 0.7583333 a
2 0.7083333 b
Comparaciones de medias para las interacciones
Para los niveles del factor A dentro del nivel B1:
A1 vs A2:
\(H_0: \mu_{A1} - \mu_{A2} = 0\)
\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)
Para los niveles del factor B dentro del nivel A1:
B1 vs B2:
\(H_0: \mu_{B1} - \mu_{B2} = 0\)
\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)
NOTA: Repetir este proceso para cada nivel de A y cada nivel de B.
Análisis de varianza para interacción de dos factores con el paquete phia
Comparación de los niveles de B dentro de cada nivel de A
multcomp.test_2factors <-function(object, respuesta, factor_name1, factor_name2, test, aov){# Función auxiliar para aplicar la prueba de comparaciones múltiples a un data frame multcomp_df <-function(df, respuesta, factor_name, test, aov){if (test =="LSD") { comp <-LSD.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[5]] } elseif (test =="HSD") { comp <-HSD.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[5]] } elseif (test =="duncan") { comp <-duncan.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[6]] } elseif (test =="SNK") { comp <-SNK.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[6]] } else {stop("Test no válido. Los tests disponibles son LSD, HSD, duncan, y SNK.") }return(comp %>%rename("y"="df[[respuesta]]") ) }# Aplicar la prueba de comparaciones múltiples a todos los data frames dentro de filters1 comp_filters1 <-lapply(object[[1]], function(df){multcomp_df(df, respuesta, factor_name2, test, aov) %>%arrange(row.names(.)) %>%rownames_to_column(var ="x") %>%relocate(x) })# Aplicar la prueba de comparaciones múltiples a todos los data frames dentro de filters2 comp_filters2 <-lapply(object[[2]], function(df){multcomp_df(df, respuesta, factor_name1, test, aov) %>% dplyr::arrange(row.names(.)) %>%rownames_to_column(var ="x") %>%relocate(x) %>% dplyr::mutate(groups =toupper(groups)) })row.names(comp_filters1) <-NULLrow.names(comp_filters2) <-NULL# Retornar una lista con las comparaciones múltiples para cada data frame result <-list() result[[as.name(substitute(factor_name1))]] <- comp_filters1 result[[as.name(substitute(factor_name2))]] <- comp_filters2return(#list(#"Comparación de los niveles del factor B dentro de cada nivel del factor A", result#[[1]],# "Comparación de los niveles del factor A dentro de cada nivel del factor B",# result[[2]]# ) )}
$Zona
$Zona$Camana
x y groups
1 1 0.7708333 a
2 2 0.6833333 b
$Zona$Castilla
x y groups
1 1 0.7458333 a
2 2 0.7333333 a
$VARIEDAD
$VARIEDAD$`1`
x y groups
1 Camana 0.7708333 A
2 Castilla 0.7458333 A
$VARIEDAD$`2`
x y groups
1 Camana 0.6833333 B
2 Castilla 0.7333333 A
# Convertir la lista a un data framedf <- result.comp %>% reshape2::melt() %>%ungroup() %>% dplyr::select(-variable) %>%relocate("Factor"="L1","Nivel"="L2", x,"y"="value", groups)df
Factor Nivel x y groups
1 Zona Camana 1 0.7708333 a
2 Zona Camana 2 0.6833333 b
3 Zona Castilla 1 0.7458333 a
4 Zona Castilla 2 0.7333333 a
5 VARIEDAD 1 Camana 0.7708333 A
6 VARIEDAD 1 Castilla 0.7458333 A
7 VARIEDAD 2 Camana 0.6833333 B
8 VARIEDAD 2 Castilla 0.7333333 A
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# # Convertir level1 y level2 en nombres simbólicos# level1 <- as.name(level1)# level2 <- as.name(level2)# # filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1, level2)) %>% dplyr::mutate(groups =paste0(groups.y,groups.x)) %>% dplyr::select(!!level1,!!level2, y, groups) #%>%# rename(!!level1 := x,# !!level2 := Nivel)# devolver una lista con los dos subconjuntosreturn(df)}
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1,level2)) %>% dplyr::mutate(y =paste0(round(y,2), " ", groups.y, groups.x)) %>% dplyr::select(!!level1,!!level2, y)# devolver una lista con los dos subconjuntosreturn(df)}
Análisis de DCA de tres vías repetido en el espacio
Definición del modelo
modelo.dca1 <-lm(Diametro_tallo_90dds ~ BLOQUE %in% Zona + DOSIS*VARIEDAD*Zona, data = data)modelo.dca2 <-lm(Diametro_tallo_90dds ~ BLOQUE %in% Zona + DOSIS+VARIEDAD*Zona, data = data)modelo.dca3 <-lm(Diametro_tallo_90dds ~ BLOQUE %in% Zona + DOSIS*Zona+VARIEDAD, data = data)modelo.dca4 <-lm(Diametro_tallo_90dds ~ BLOQUE %in% Zona + DOSIS*VARIEDAD+Zona, data = data)modelo.dca5 <-lm(Diametro_tallo_90dds ~ BLOQUE %in% Zona + DOSIS+VARIEDAD+Zona, data = data)
Nota
La expresión “BLOQUE/ZONA” es para considerar la interacción como efecto fijo.
La expresión “BLOQUE:ZONA” es para considerar la interacción como efecto aleatorio.
La expresión “BLOQUE %in% ZONA” permite que se genere una estructura del análisis de varianza donde se anide el efecto de la variedad dentro de cada bloque
Durbin-Watson test
data: modelo.dca
DW = 2.2589, p-value = 0.6066
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
shapiro.test(rstudent(modelo.dca))
Shapiro-Wilk normality test
data: rstudent(modelo.dca)
W = 0.96809, p-value = 0.2132
ad.test(rstudent(modelo.dca))
Anderson-Darling normality test
data: rstudent(modelo.dca)
A = 0.51434, p-value = 0.1832
lillie.test(rstudent(modelo.dca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dca)
D = 0.092111, p-value = 0.3907
Asymptotic one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dca)
D = 0.09263, p-value = 0.8047
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dca))
Cramer-von Mises normality test
data: rstudent(modelo.dca)
W = 0.079806, p-value = 0.2033
pearson.test(rstudent(modelo.dca))
Pearson chi-square normality test
data: rstudent(modelo.dca)
P = 10.333, p-value = 0.1705
sf.test(rstudent(modelo.dca))
Shapiro-Francia normality test
data: rstudent(modelo.dca)
W = 0.96394, p-value = 0.1295
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.7436697, Df = 1, p = 0.38849
bptest(modelo.dca)
studentized Breusch-Pagan test
data: modelo.dca
BP = 11.834, df = 9, p-value = 0.2228
bptest(modelo.dca, studentize = F)
Breusch-Pagan test
data: modelo.dca
BP = 14.449, df = 9, p-value = 0.1072
olsrr::ols_test_breusch_pagan(modelo.dca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
------------------------------------------------
Response : Diametro_tallo_90dds
Variables: fitted values of Diametro_tallo_90dds
Test Summary
----------------------------
DF = 1
Chi2 = 0.7436697
Prob > Chi2 = 0.3884879
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.
Estadísticas globales
modelo.dca %>%gvlma()
Call:
lm(formula = Diametro_tallo_90dds ~ BLOQUE %in% Zona + DOSIS +
VARIEDAD + Zona, data = data)
Coefficients:
(Intercept) DOSIS1 DOSIS2
7.375e-01 1.667e-02 1.667e-02
DOSIS3 VARIEDAD2 ZonaCastilla
1.667e-02 -5.000e-02 -6.250e-03
BLOQUEII:ZonaCamana BLOQUEIII:ZonaCamana BLOQUEII:ZonaCastilla
-1.250e-02 1.875e-02 6.250e-02
BLOQUEIII:ZonaCastilla
-2.602e-18
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = .)
Value p-value Decision
Global Stat 1.804587 0.7716 Assumptions acceptable.
Skewness 0.181520 0.6701 Assumptions acceptable.
Kurtosis 0.390464 0.5321 Assumptions acceptable.
Link Function 1.228629 0.2677 Assumptions acceptable.
Heteroscedasticity 0.003974 0.9497 Assumptions acceptable.
\(Y_{ijkl}\) = Valor observado de la variable respuesta.
\(\hat{Y}_{ijkl}\) = Valor ajustado de la variable respuesta.
\(\mu\) = Promedio observado de la variable respuesta.
\(\tau_{i}\) = Efecto del i-ésimo nivel del factor Zona.
\(\text{Error}(\tau\text{rep})_{i(k)}\) = Efecto del i-ésimo nivel de Zona en el k-ésimo nivel de repetición.
\(\beta_{j}\) = Efecto del j-ésimo nivel del factor tratamiento (Variedad).
\(\gamma_{k}\) = Efecto del k-ésimo nivel del factor tratamiento (Dosis).
\(\epsilon_{ijkl}\) = Residuo observado del modelo.
Pruebas de hipótesis
Para el factor A (Zona):
\(H_0: \tau_{A1} = \tau_{A2} = 0\)
\(H_1: \text{En al menos un nivel del factor A el } \tau \text{ es diferente a los demás.}\)
\(H_1: \tau_i \neq 0\text{; en al menos un nivel del factor A.}\)
Para el factor B (Variedad):
\(H_0: \beta_{B1} = \beta_{B2} = 0\)
\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)
\(H_1: \beta_j \neq 0\text{; en al menos un nivel del factor B.}\)
Para el factor C (Dosis):
\(H_0: \gamma_{C1} = \gamma_{C2} = 0\)
\(H_1: \text{En al menos un nivel del factor C el } \beta \text{ es diferente a los demás.}\)
\(H_1: \gamma_k \neq 0\text{; en al menos un nivel del factor C.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
anova(modelo.dca, test ="F")
Analysis of Variance Table
Response: Diametro_tallo_90dds
Df Sum Sq Mean Sq F value Pr(>F)
DOSIS 3 0.002500 0.0008333 0.1306 0.94131
VARIEDAD 1 0.030000 0.0300000 4.7010 0.03647 *
Zona 1 0.001875 0.0018750 0.2938 0.59095
BLOQUE:Zona 4 0.024792 0.0061979 0.9712 0.43458
Residuals 38 0.242500 0.0063816
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
GAD::estimates(modelo.dca)
$tm
BLOQUE Zona DOSIS VARIEDAD n
DOSIS 3 2 0 2 4.8
VARIEDAD 3 2 4 0 4.8
Zona 3 0 4 2 4.8
BLOQUE:Zona 1 1 4 2 4.8
Res 1 1 1 1 1.0
$mse
Mean square estimates
DOSIS "Res + DOSIS"
VARIEDAD "Res + VARIEDAD"
Zona "Res + BLOQUE:Zona + Zona"
BLOQUE:Zona "Res + BLOQUE:Zona"
Residual "Res"
$f.versus
F-ratio versus
DOSIS "Residual"
VARIEDAD "Residual"
Zona "BLOQUE:Zona"
BLOQUE:Zona "Residual"
GAD::gad(modelo.dca)
Analysis of Variance Table
Response: Diametro_tallo_90dds
Df Sum Sq Mean Sq F value Pr(>F)
DOSIS 3 0.002500 0.0008333 0.1306 0.94131
VARIEDAD 1 0.030000 0.0300000 4.7010 0.03647 *
Zona 1 0.001875 0.0018750 0.3025 0.61157
BLOQUE:Zona 4 0.024792 0.0061979 0.9712 0.43458
Residual 38 0.242500 0.0063816
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Importante
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(Zona/rep)” en el caso de efectos aleatorios o “Error(Zona:rep)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos.
aov(Diametro_tallo_90dds ~ Zona + VARIEDAD + DOSIS +Error(BLOQUE %in% Zona), data = data) -> aov.dcasummary(aov.dca)
Error: BLOQUE:Zona
Df Sum Sq Mean Sq F value Pr(>F)
Zona 1 0.001875 0.001875 0.303 0.612
Residuals 4 0.024792 0.006198
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
VARIEDAD 1 0.0300 0.030000 4.701 0.0365 *
DOSIS 3 0.0025 0.000833 0.131 0.9413
Residuals 38 0.2425 0.006382
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
broom::tidy(aov.dca)
# A tibble: 5 × 7
stratum term df sumsq meansq statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 BLOQUE:Zona Zona 1 0.00187 0.00187 0.303 0.612
2 BLOQUE:Zona Residuals 4 0.0248 0.00620 NA NA
3 Within VARIEDAD 1 0.0300 0.0300 4.70 0.0365
4 Within DOSIS 3 0.00250 0.000833 0.131 0.941
5 Within Residuals 38 0.243 0.00638 NA NA
Valor de la tabla de F para el factor Zona con una significancia de 0.05.
qf(0.95, 1, 4)
[1] 7.708647
Valor de la tabla de F para el factor Variedad con una significancia de 0.05.
qf(0.95, 1, 38)
[1] 4.098172
Valor de la tabla de F para el factor Dosis con una significancia de 0.05.
qf(0.95, 3, 38)
[1] 2.851741
Conclusión.
Con respecto al Factor Zona: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor Zona tienen un efecto sobre el diámetro de tallo estadísticamente similar a 0.
Con respecto al Factor Variedad: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor Variedad tienen un efecto sobre el diámetro de tallo estadísticamente diferente a 0.
Con respecto al Factor Dosis: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor Dosis tienen un efecto sobre el diámetro de tallo estadísticamente similar a 0.
agricolae::cv.model(modelo.dca)
[1] 10.89338
Comparaciones de medias para los efectos principales del Factor Zona
data %>%with(LSD.test( Diametro_tallo_90dds, # Cambiar según nombre de variable respuesta Zona, # Cambiar según nombre de variable independienteDFerror =get_df_ea(aov.dca), MSerror =get_mse_ea(aov.dca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: Diametro_tallo_90dds ~ Zona
LSD t Test for Diametro_tallo_90dds
Mean Square Error: 0.006197917
Zona, means and individual ( 95 %) CI
Diametro_tallo_90dds std r LCL UCL Min Max
Camana 0.7270833 0.09205288 24 0.6824657 0.7717009 0.5 0.90
Castilla 0.7395833 0.06753287 24 0.6949657 0.7842009 0.6 0.85
Alpha: 0.05 ; DF Error: 4
Critical Value of t: 2.776445
least Significant Difference: 0.06309883
Treatments with the same letter are not significantly different.
Diametro_tallo_90dds groups
Castilla 0.7395833 a
Camana 0.7270833 a
Comparaciones de medias para los efectos principales del Factor Variedad
data %>%with(LSD.test( Diametro_tallo_90dds, # Cambiar según nombre de variable respuesta VARIEDAD, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dca), MSerror = dvmisc::get_mse(modelo.dca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: Diametro_tallo_90dds ~ VARIEDAD
LSD t Test for Diametro_tallo_90dds
Mean Square Error: 0.006381579
VARIEDAD, means and individual ( 95 %) CI
Diametro_tallo_90dds std r LCL UCL Min Max
1 0.7583333 0.07322786 24 0.7253227 0.7913439 0.6 0.90
2 0.7083333 0.08030738 24 0.6753227 0.7413439 0.5 0.85
Alpha: 0.05 ; DF Error: 38
Critical Value of t: 2.024394
least Significant Difference: 0.04668405
Treatments with the same letter are not significantly different.
Diametro_tallo_90dds groups
1 0.7583333 a
2 0.7083333 b
Comparaciones de medias para los efectos principales del Factor Dosis
data %>%with(LSD.test( Diametro_tallo_90dds, # Cambiar según nombre de variable respuesta DOSIS, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dca), MSerror = dvmisc::get_mse(modelo.dca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: Diametro_tallo_90dds ~ DOSIS
LSD t Test for Diametro_tallo_90dds
Mean Square Error: 0.006381579
DOSIS, means and individual ( 95 %) CI
Diametro_tallo_90dds std r LCL UCL Min Max
0 0.7208333 0.07216878 12 0.6741493 0.7675174 0.6 0.8
1 0.7375000 0.06440285 12 0.6908159 0.7841841 0.6 0.8
2 0.7375000 0.06784005 12 0.6908159 0.7841841 0.6 0.8
3 0.7375000 0.11505927 12 0.6908159 0.7841841 0.5 0.9
Alpha: 0.05 ; DF Error: 38
Critical Value of t: 2.024394
least Significant Difference: 0.06602122
Treatments with the same letter are not significantly different.
Diametro_tallo_90dds groups
1 0.7375000 a
2 0.7375000 a
3 0.7375000 a
0 0.7208333 a
Análisis de DBCA de tres vías repetido en el espacio
Definición del modelo
modelo.dbca1 <-lm(Diametro_tallo_90dds ~ BLOQUE + BLOQUE %in% Zona + DOSIS*VARIEDAD*Zona, data = data)modelo.dbca2 <-lm(Diametro_tallo_90dds ~ BLOQUE + BLOQUE %in% Zona + DOSIS+VARIEDAD*Zona, data = data)modelo.dbca3 <-lm(Diametro_tallo_90dds ~ BLOQUE + BLOQUE %in% Zona + DOSIS*Zona+VARIEDAD, data = data)modelo.dbca4 <-lm(Diametro_tallo_90dds ~ BLOQUE + BLOQUE %in% Zona + DOSIS*VARIEDAD+Zona, data = data)modelo.dbca5 <-lm(Diametro_tallo_90dds ~ BLOQUE + BLOQUE %in% Zona + DOSIS+VARIEDAD+Zona, data = data)
Nota
La expresión “BLOQUE/ZONA” es para considerar la interacción como efecto fijo.
La expresión “BLOQUE:ZONA” es para considerar la interacción como efecto aleatorio.
La expresión “BLOQUE %in% ZONA” permite que se genere una estructura del análisis de varianza donde se anide el efecto de la variedad dentro de cada bloque
Durbin-Watson test
data: modelo.dbca
DW = 2.8708, p-value = 0.2701
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
shapiro.test(rstudent(modelo.dbca))
Shapiro-Wilk normality test
data: rstudent(modelo.dbca)
W = 0.96824, p-value = 0.2162
ad.test(rstudent(modelo.dbca))
Anderson-Darling normality test
data: rstudent(modelo.dbca)
A = 0.5146, p-value = 0.183
lillie.test(rstudent(modelo.dbca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dbca)
D = 0.14251, p-value = 0.01599
Asymptotic one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dbca)
D = 0.14768, p-value = 0.246
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dbca))
Cramer-von Mises normality test
data: rstudent(modelo.dbca)
W = 0.086058, p-value = 0.1685
pearson.test(rstudent(modelo.dbca))
Pearson chi-square normality test
data: rstudent(modelo.dbca)
P = 11.167, p-value = 0.1315
sf.test(rstudent(modelo.dbca))
Shapiro-Francia normality test
data: rstudent(modelo.dbca)
W = 0.97236, p-value = 0.2651
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dbca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.10062, Df = 1, p = 0.75109
bptest(modelo.dbca)
studentized Breusch-Pagan test
data: modelo.dbca
BP = 24.394, df = 19, p-value = 0.1815
bptest(modelo.dbca, studentize = F)
Breusch-Pagan test
data: modelo.dbca
BP = 15.27, df = 19, p-value = 0.7053
olsrr::ols_test_breusch_pagan(modelo.dbca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
------------------------------------------------
Response : Diametro_tallo_90dds
Variables: fitted values of Diametro_tallo_90dds
Test Summary
----------------------------
DF = 1
Chi2 = 0.10062
Prob > Chi2 = 0.7510868
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.
Estadísticas globales
modelo.dbca %>%gvlma()
Call:
lm(formula = Diametro_tallo_90dds ~ BLOQUE + BLOQUE %in% Zona +
DOSIS * VARIEDAD * Zona, data = data)
Coefficients:
(Intercept) BLOQUEII
7.646e-01 -1.250e-02
BLOQUEIII DOSIS1
1.875e-02 -1.381e-15
DOSIS2 DOSIS3
-5.000e-02 6.667e-02
VARIEDAD2 ZonaCastilla
-5.000e-02 -6.875e-02
BLOQUEII:ZonaCastilla BLOQUEIII:ZonaCastilla
7.500e-02 -1.875e-02
DOSIS1:VARIEDAD2 DOSIS2:VARIEDAD2
-5.000e-02 6.667e-02
DOSIS3:VARIEDAD2 ZonaCastilla:DOSIS1
-1.667e-01 5.000e-02
ZonaCastilla:DOSIS2 ZonaCastilla:DOSIS3
1.167e-01 -6.667e-02
ZonaCastilla:VARIEDAD2 ZonaCastilla:DOSIS1:VARIEDAD2
1.667e-02 6.667e-02
ZonaCastilla:DOSIS2:VARIEDAD2 ZonaCastilla:DOSIS3:VARIEDAD2
-1.000e-01 2.667e-01
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = .)
Value p-value Decision
Global Stat 10.2068 0.037084 Assumptions NOT satisfied!
Skewness 0.7767 0.378153 Assumptions acceptable.
Kurtosis 1.1192 0.290098 Assumptions acceptable.
Link Function 8.1303 0.004353 Assumptions NOT satisfied!
Heteroscedasticity 0.1807 0.670787 Assumptions acceptable.
\(Y_{ijkl}\) = Valor observado de la variable respuesta.
\(\hat{Y}_{ijkl}\) = Valor ajustado de la variable respuesta.
\(\mu\) = Promedio observado de la variable respuesta.
\(\tau_{i}\) = Efecto del i-ésimo nivel del factor Zona.
\(\text{Error}(\tau\text{Bloque})_{i(l)}\) = Efecto del i-ésimo nivel de Zona en el l-ésimo nivel de Bloque.
\(\beta_{j}\) = Efecto del j-ésimo nivel del factor tratamiento (Variedad).
\(\gamma_{k}\) = Efecto del k-ésimo nivel del factor tratamiento (Dosis).
\(\delta_{l}\) = Efecto del l-ésimo nivel del factor Bloque.
\(\epsilon_{ijkl}\) = Residuo observado del modelo.
Pruebas de hipótesis
Para el factor A (Zona):
\(H_0: \tau_{A1} = \tau_{A2} = 0\)
\(H_1: \text{En al menos un nivel del factor A el } \tau \text{ es diferente a los demás.}\)
\(H_1: \tau_i \neq 0\text{; en al menos un nivel del factor A.}\)
Para el factor B (Variedad):
\(H_0: \beta_{B1} = \beta_{B2} = 0\)
\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)
\(H_1: \beta_j \neq 0\text{; en al menos un nivel del factor B.}\)
Para el factor C (Dosis):
\(H_0: \gamma_{C1} = \gamma_{C2} = 0\)
\(H_1: \text{En al menos un nivel del factor C el } \beta \text{ es diferente a los demás.}\)
\(H_1: \gamma_k \neq 0\text{; en al menos un nivel del factor C.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(Zona/Bloque)” en el caso de efectos aleatorios o “Error(Zona:Bloque)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos.
aov(Diametro_tallo_90dds ~ BLOQUE + Zona * VARIEDAD * DOSIS +Error(BLOQUE %in% Zona), data = data) -> aov.dbcasummary(aov.dbca)
Error: BLOQUE:Zona
Df Sum Sq Mean Sq F value Pr(>F)
BLOQUE 2 0.005104 0.002552 0.259 0.794
Zona 1 0.001875 0.001875 0.190 0.705
Residuals 2 0.019688 0.009844
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
VARIEDAD 1 0.03000 0.030000 5.412 0.0275 *
DOSIS 3 0.00250 0.000833 0.150 0.9286
Zona:VARIEDAD 1 0.01688 0.016875 3.044 0.0920 .
Zona:DOSIS 3 0.01229 0.004097 0.739 0.5376
VARIEDAD:DOSIS 3 0.00417 0.001389 0.251 0.8602
Zona:VARIEDAD:DOSIS 3 0.05396 0.017986 3.245 0.0368 *
Residuals 28 0.15521 0.005543
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
broom::tidy(aov.dbca)
# A tibble: 10 × 7
stratum term df sumsq meansq statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 BLOQUE:Zona BLOQUE 2 0.00510 0.00255 0.259 0.794
2 BLOQUE:Zona Zona 1 0.00187 0.00187 0.190 0.705
3 BLOQUE:Zona Residuals 2 0.0197 0.00984 NA NA
4 Within VARIEDAD 1 0.0300 0.0300 5.41 0.0275
5 Within DOSIS 3 0.00250 0.000833 0.150 0.929
6 Within Zona:VARIEDAD 1 0.0169 0.0169 3.04 0.0920
7 Within Zona:DOSIS 3 0.0123 0.00410 0.739 0.538
8 Within VARIEDAD:DOSIS 3 0.00417 0.00139 0.251 0.860
9 Within Zona:VARIEDAD:DOSIS 3 0.0540 0.0180 3.24 0.0368
10 Within Residuals 28 0.155 0.00554 NA NA
Valor de la tabla de F para el factor Zona con una significancia de 0.05.
qf(0.95, 1, 2)
[1] 18.51282
Valor de la tabla de F para el factor Variedad con una significancia de 0.05.
qf(0.95, 1, 28)
[1] 4.195972
Valor de la tabla de F para el factor Dosis con una significancia de 0.05.
qf(0.95, 3, 28)
[1] 2.946685
Conclusión.
Con respecto al Factor Zona: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor Zona tienen un efecto sobre el diámetro de tallo estadísticamente similar a 0.
Con respecto al Factor Variedad: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor Variedad tienen un efecto sobre el diámetro de tallo estadísticamente similar a 0.
Con respecto al Factor Dosis: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor Dosis tienen un efecto sobre el diámetro de tallo estadísticamente similar a 0.
Con respecto a la interacción entre Zona/Variedad/Dosis: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, existe al menos un nivel de interacción estadísticamente antagonista y sinergista con respecto al diámetro de tallo.
agricolae::cv.model(modelo.dbca)
[1] 10.1526
Comparaciones de medias para los efectos principales del Factor Zona
data %>%with(LSD.test( Diametro_tallo_90dds, # Cambiar según nombre de variable respuesta Zona, # Cambiar según nombre de variable independienteDFerror =get_df_ea(aov.dbca), MSerror =get_mse_ea(aov.dbca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: Diametro_tallo_90dds ~ Zona
LSD t Test for Diametro_tallo_90dds
Mean Square Error: 0.00984375
Zona, means and individual ( 95 %) CI
Diametro_tallo_90dds std r LCL UCL Min Max
Camana 0.7270833 0.09205288 24 0.6399447 0.814222 0.5 0.90
Castilla 0.7395833 0.06753287 24 0.6524447 0.826722 0.6 0.85
Alpha: 0.05 ; DF Error: 2
Critical Value of t: 4.302653
least Significant Difference: 0.1232327
Treatments with the same letter are not significantly different.
Diametro_tallo_90dds groups
Castilla 0.7395833 a
Camana 0.7270833 a
Comparaciones de medias para los efectos principales del Factor Variedad
data %>%with(LSD.test( Diametro_tallo_90dds, # Cambiar según nombre de variable respuesta VARIEDAD, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dbca), MSerror = dvmisc::get_mse(modelo.dbca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: Diametro_tallo_90dds ~ VARIEDAD
LSD t Test for Diametro_tallo_90dds
Mean Square Error: 0.005543155
VARIEDAD, means and individual ( 95 %) CI
Diametro_tallo_90dds std r LCL UCL Min Max
1 0.7583333 0.07322786 24 0.7272026 0.7894641 0.6 0.90
2 0.7083333 0.08030738 24 0.6772026 0.7394641 0.5 0.85
Alpha: 0.05 ; DF Error: 28
Critical Value of t: 2.048407
least Significant Difference: 0.04402549
Treatments with the same letter are not significantly different.
Diametro_tallo_90dds groups
1 0.7583333 a
2 0.7083333 b
Comparaciones de medias para las interacciones
Para cada Zona, se debe aplicar lo siguiente:
Para los niveles del factor B dentro del nivel C1:
B1 vs B2:
\(H_0: \mu_{B1} - \mu_{B2} = 0\)
\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)
…
B3 vs B4:
\(H_0: \mu_{B3} - \mu_{B4} = 0\)
\(H_1: \mu_{B4} - \mu_{B4} \neq 0\)
Para los niveles del factor C dentro del nivel B1:
C1 vs C2:
\(H_0: \mu_{C1} - \mu_{C2} = 0\)
\(H_1: \mu_{C1} - \mu_{C2} \neq 0\)
NOTA: Repetir este proceso para cada nivel de B y cada nivel de C.
Análisis de varianza para interacción de tres factores con el paquete emmeans
Comparación de la interacción de los niveles de B y los niveles de C dentro de cada nivel de A
emmeans::joint_tests(modelo.dbca, by ="Zona")
Zona = Camana:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 0.357 0.7029
DOSIS 3 28 0.119 0.9482
VARIEDAD 1 28 8.287 0.0076
DOSIS:VARIEDAD 3 28 2.625 0.0701
Zona = Castilla:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 1.879 0.1715
DOSIS 3 28 0.770 0.5203
VARIEDAD 1 28 0.169 0.6840
DOSIS:VARIEDAD 3 28 0.871 0.4679
emmeans::joint_tests(modelo.dbca, by =c("Zona","VARIEDAD"))
Zona = Camana, VARIEDAD = 1:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 0.357 0.7029
DOSIS 3 28 1.240 0.3138
Zona = Castilla, VARIEDAD = 1:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 1.879 0.1715
DOSIS 3 28 0.639 0.5963
Zona = Camana, VARIEDAD = 2:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 0.357 0.7029
DOSIS 3 28 1.503 0.2353
Zona = Castilla, VARIEDAD = 2:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 1.879 0.1715
DOSIS 3 28 1.002 0.4064
emmeans::joint_tests(modelo.dbca, by =c("Zona","DOSIS"))
Zona = Camana, DOSIS = 0:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 0.357 0.7029
VARIEDAD 1 28 0.677 0.4177
Zona = Castilla, DOSIS = 0:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 1.879 0.1715
VARIEDAD 1 28 0.301 0.5878
Zona = Camana, DOSIS = 1:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 0.357 0.7029
VARIEDAD 1 28 2.706 0.1112
Zona = Castilla, DOSIS = 1:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 1.879 0.1715
VARIEDAD 1 28 0.075 0.7860
Zona = Camana, DOSIS = 2:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 0.357 0.7029
VARIEDAD 1 28 0.075 0.7860
Zona = Castilla, DOSIS = 2:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 1.879 0.1715
VARIEDAD 1 28 1.203 0.2821
Zona = Camana, DOSIS = 3:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 0.357 0.7029
VARIEDAD 1 28 12.703 0.0013
Zona = Castilla, DOSIS = 3:
model term df1 df2 F.ratio p.value
BLOQUE 2 28 1.879 0.1715
VARIEDAD 1 28 1.203 0.2821
DOSIS = 0, Zona = Camana:
VARIEDAD emmean SE df lower.CL upper.CL .group
2 0.717 0.043 28 0.629 0.805 A
1 0.767 0.043 28 0.679 0.855 A
DOSIS = 1, Zona = Camana:
VARIEDAD emmean SE df lower.CL upper.CL .group
2 0.667 0.043 28 0.579 0.755 A
1 0.767 0.043 28 0.679 0.855 A
DOSIS = 2, Zona = Camana:
VARIEDAD emmean SE df lower.CL upper.CL .group
1 0.717 0.043 28 0.629 0.805 A
2 0.733 0.043 28 0.645 0.821 A
DOSIS = 3, Zona = Camana:
VARIEDAD emmean SE df lower.CL upper.CL .group
2 0.617 0.043 28 0.529 0.705 A
1 0.833 0.043 28 0.745 0.921 B
DOSIS = 0, Zona = Castilla:
VARIEDAD emmean SE df lower.CL upper.CL .group
2 0.683 0.043 28 0.595 0.771 A
1 0.717 0.043 28 0.629 0.805 A
DOSIS = 1, Zona = Castilla:
VARIEDAD emmean SE df lower.CL upper.CL .group
2 0.750 0.043 28 0.662 0.838 A
1 0.767 0.043 28 0.679 0.855 A
DOSIS = 2, Zona = Castilla:
VARIEDAD emmean SE df lower.CL upper.CL .group
2 0.717 0.043 28 0.629 0.805 A
1 0.783 0.043 28 0.695 0.871 A
DOSIS = 3, Zona = Castilla:
VARIEDAD emmean SE df lower.CL upper.CL .group
1 0.717 0.043 28 0.629 0.805 A
2 0.783 0.043 28 0.695 0.871 A
Results are averaged over the levels of: BLOQUE
Confidence level used: 0.95
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
VARIEDAD = 1, Zona = Camana:
DOSIS emmean SE df lower.CL upper.CL .group
2 0.717 0.043 28 0.629 0.805 A
1 0.767 0.043 28 0.679 0.855 A
0 0.767 0.043 28 0.679 0.855 A
3 0.833 0.043 28 0.745 0.921 A
VARIEDAD = 2, Zona = Camana:
DOSIS emmean SE df lower.CL upper.CL .group
3 0.617 0.043 28 0.529 0.705 A
1 0.667 0.043 28 0.579 0.755 A
0 0.717 0.043 28 0.629 0.805 A
2 0.733 0.043 28 0.645 0.821 A
VARIEDAD = 1, Zona = Castilla:
DOSIS emmean SE df lower.CL upper.CL .group
0 0.717 0.043 28 0.629 0.805 A
3 0.717 0.043 28 0.629 0.805 A
1 0.767 0.043 28 0.679 0.855 A
2 0.783 0.043 28 0.695 0.871 A
VARIEDAD = 2, Zona = Castilla:
DOSIS emmean SE df lower.CL upper.CL .group
0 0.683 0.043 28 0.595 0.771 A
2 0.717 0.043 28 0.629 0.805 A
1 0.750 0.043 28 0.662 0.838 A
3 0.783 0.043 28 0.695 0.871 A
Results are averaged over the levels of: BLOQUE
Confidence level used: 0.95
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
niveles_zona <-unique(data$Zona)# Aplicar la función filter_by_2factor_level a cada nivel de Zonadatos_filtrados <- purrr::map( niveles_zona, ~filter_by_2factor_level(data =filter(data, Zona == .x),factor_name1 = VARIEDAD,factor_name2 = DOSIS))names(datos_filtrados) <- niveles_zona
$Camana
$Camana$VARIEDAD
$Camana$VARIEDAD$`1`
x y groups
1 0 0.7666667 a
2 1 0.7666667 a
3 2 0.7166667 a
4 3 0.8333333 a
$Camana$VARIEDAD$`2`
x y groups
1 0 0.7166667 a
2 1 0.6666667 a
3 2 0.7333333 a
4 3 0.6166667 a
$Camana$DOSIS
$Camana$DOSIS$`0`
x y groups
1 1 0.7666667 A
2 2 0.7166667 A
$Camana$DOSIS$`1`
x y groups
1 1 0.7666667 A
2 2 0.6666667 A
$Camana$DOSIS$`2`
x y groups
1 1 0.7166667 A
2 2 0.7333333 A
$Camana$DOSIS$`3`
x y groups
1 1 0.8333333 A
2 2 0.6166667 B
$Castilla
$Castilla$VARIEDAD
$Castilla$VARIEDAD$`1`
x y groups
1 0 0.7166667 a
2 1 0.7666667 a
3 2 0.7833333 a
4 3 0.7166667 a
$Castilla$VARIEDAD$`2`
x y groups
1 0 0.6833333 a
2 1 0.7500000 a
3 2 0.7166667 a
4 3 0.7833333 a
$Castilla$DOSIS
$Castilla$DOSIS$`0`
x y groups
1 1 0.7166667 A
2 2 0.6833333 A
$Castilla$DOSIS$`1`
x y groups
1 1 0.7666667 A
2 2 0.7500000 A
$Castilla$DOSIS$`2`
x y groups
1 1 0.7833333 A
2 2 0.7166667 A
$Castilla$DOSIS$`3`
x y groups
1 1 0.7166667 A
2 2 0.7833333 A
# Convertir la lista a un data framedf <- result.comp %>% reshape2::melt() %>%ungroup() %>% dplyr::select(-variable) %>%relocate("Zona"="L1","Factor"="L2","Nivel"="L3", x,"y"="value", groups)df
Zona Factor Nivel x y groups
1 Camana VARIEDAD 1 0 0.7666667 a
2 Camana VARIEDAD 1 1 0.7666667 a
3 Camana VARIEDAD 1 2 0.7166667 a
4 Camana VARIEDAD 1 3 0.8333333 a
5 Camana VARIEDAD 2 0 0.7166667 a
6 Camana VARIEDAD 2 1 0.6666667 a
7 Camana VARIEDAD 2 2 0.7333333 a
8 Camana VARIEDAD 2 3 0.6166667 a
9 Camana DOSIS 0 1 0.7666667 A
10 Camana DOSIS 0 2 0.7166667 A
11 Camana DOSIS 1 1 0.7666667 A
12 Camana DOSIS 1 2 0.6666667 A
13 Camana DOSIS 2 1 0.7166667 A
14 Camana DOSIS 2 2 0.7333333 A
15 Camana DOSIS 3 1 0.8333333 A
16 Camana DOSIS 3 2 0.6166667 B
17 Castilla VARIEDAD 1 0 0.7166667 a
18 Castilla VARIEDAD 1 1 0.7666667 a
19 Castilla VARIEDAD 1 2 0.7833333 a
20 Castilla VARIEDAD 1 3 0.7166667 a
21 Castilla VARIEDAD 2 0 0.6833333 a
22 Castilla VARIEDAD 2 1 0.7500000 a
23 Castilla VARIEDAD 2 2 0.7166667 a
24 Castilla VARIEDAD 2 3 0.7833333 a
25 Castilla DOSIS 0 1 0.7166667 A
26 Castilla DOSIS 0 2 0.6833333 A
27 Castilla DOSIS 1 1 0.7666667 A
28 Castilla DOSIS 1 2 0.7500000 A
29 Castilla DOSIS 2 1 0.7833333 A
30 Castilla DOSIS 2 2 0.7166667 A
31 Castilla DOSIS 3 1 0.7166667 A
32 Castilla DOSIS 3 2 0.7833333 A
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# # Convertir level1 y level2 en nombres simbólicos# level1 <- as.name(level1)# level2 <- as.name(level2)# # filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1, level2)) %>% dplyr::mutate(groups =paste0(groups.y,groups.x)) %>% dplyr::select(!!level1,!!level2, y, groups) #%>%# rename(!!level1 := x,# !!level2 := Nivel)# devolver una lista con los dos subconjuntosreturn(df)}
zonas <-unique(df$Zona)resultados <- zonas %>%map(~create_report(filter(df, Zona == .x), "VARIEDAD", "DOSIS"))names(resultados) <- zonasresultados
$Camana
VARIEDAD DOSIS y groups
1 1 0 0.7666667 Aa
2 1 1 0.7666667 Aa
3 1 2 0.7166667 Aa
4 1 3 0.8333333 Aa
5 2 0 0.7166667 Aa
6 2 1 0.6666667 Aa
7 2 2 0.7333333 Aa
8 2 3 0.6166667 Ba
$Castilla
VARIEDAD DOSIS y groups
1 1 0 0.7166667 Aa
2 1 1 0.7666667 Aa
3 1 2 0.7833333 Aa
4 1 3 0.7166667 Aa
5 2 0 0.6833333 Aa
6 2 1 0.7500000 Aa
7 2 2 0.7166667 Aa
8 2 3 0.7833333 Aa
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1,level2)) %>% dplyr::mutate(y =paste0(round(y,2), " ", groups.y, groups.x)) %>% dplyr::select(!!level1,!!level2, y)# devolver una lista con los dos subconjuntosreturn(df)}
apply_create_report <-function(df, level1, level2, ExpRep) { zonas <-unique(df[[ExpRep]]) resultados <- purrr::map(zonas, ~create_report(filter(df, .data[[ExpRep]] == .x), level1, level2))# Asignar los nombres de las zonas a los resultadosnames(resultados) <- zonas# Combinar los resultados en un data frame final df_final <-bind_rows(resultados, .id = ExpRep)return(df_final)}
Durbin-Watson test
data: modelo.dbca
DW = 2.5533, p-value = 0.3762
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
shapiro.test(rstudent(modelo.dbca))
Shapiro-Wilk normality test
data: rstudent(modelo.dbca)
W = 0.97797, p-value = 0.497
ad.test(rstudent(modelo.dbca))
Anderson-Darling normality test
data: rstudent(modelo.dbca)
A = 0.44887, p-value = 0.2665
lillie.test(rstudent(modelo.dbca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dbca)
D = 0.10527, p-value = 0.2024
Exact one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dbca)
D = 0.098701, p-value = 0.7009
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dbca))
Cramer-von Mises normality test
data: rstudent(modelo.dbca)
W = 0.07727, p-value = 0.2195
pearson.test(rstudent(modelo.dbca))
Pearson chi-square normality test
data: rstudent(modelo.dbca)
P = 6.1667, p-value = 0.5204
sf.test(rstudent(modelo.dbca))
Shapiro-Francia normality test
data: rstudent(modelo.dbca)
W = 0.97561, p-value = 0.3478
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dbca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 11.05181, Df = 1, p = 0.00088601
bptest(modelo.dbca)
studentized Breusch-Pagan test
data: modelo.dbca
BP = 22.771, df = 17, p-value = 0.1568
bptest(modelo.dbca, studentize = F)
Breusch-Pagan test
data: modelo.dbca
BP = 24.273, df = 17, p-value = 0.1122
olsrr::ols_test_breusch_pagan(modelo.dbca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
--------------------------------
Response : rdto
Variables: fitted values of rdto
Test Summary
-------------------------------
DF = 1
Chi2 = 11.05181
Prob > Chi2 = 0.0008860059
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.
Estadísticas globales
modelo.dbca %>%gvlma()
Call:
lm(formula = rdto ~ lugar * linaje + bloque %in% lugar + bloque,
data = data_D)
Coefficients:
(Intercept) lugarPisco linaje2
11.0833 1.9167 -0.1667
linaje3 linaje4 bloqueII
-1.5000 -3.6667 1.7500
bloqueIII bloqueIV bloqueV
0.7500 1.0000 2.5000
bloqueVI lugarPisco:linaje2 lugarPisco:linaje3
2.5000 2.3333 2.5000
lugarPisco:linaje4 lugarPisco:bloqueII lugarPisco:bloqueIII
3.5000 1.2500 -1.0000
lugarPisco:bloqueIV lugarPisco:bloqueV lugarPisco:bloqueVI
4.0000 -0.5000 -1.2500
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = .)
Value p-value Decision
Global Stat 7.6911887 0.103568 Assumptions acceptable.
Skewness 0.0622679 0.802947 Assumptions acceptable.
Kurtosis 0.0347871 0.852043 Assumptions acceptable.
Link Function 0.0003691 0.984673 Assumptions acceptable.
Heteroscedasticity 7.5937648 0.005857 Assumptions NOT satisfied!
\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)
\(H_1: \beta_j \neq 0\text{; en al menos un nivel del factor B.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(Zona/Bloque)” en el caso de efectos aleatorios o “Error(Zona:Bloque)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos.
aov(rdto ~ bloque + lugar * linaje +Error(bloque %in% lugar) + bloque, data = data_D) -> aov.dbcasummary(aov.dbca)
Error: bloque:lugar
Df Sum Sq Mean Sq F value Pr(>F)
bloque 5 59.50 11.90 1.539 0.32389
lugar 1 234.08 234.08 30.269 0.00271 **
Residuals 5 38.67 7.73
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
linaje 3 52.75 17.583 1.71 0.186
lugar:linaje 3 19.75 6.583 0.64 0.595
Residuals 30 308.50 10.283
broom::tidy(aov.dbca)
# A tibble: 6 × 7
stratum term df sumsq meansq statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 bloque:lugar bloque 5 59.5 11.9 1.54 0.324
2 bloque:lugar lugar 1 234. 234. 30.3 0.00271
3 bloque:lugar Residuals 5 38.7 7.73 NA NA
4 Within linaje 3 52.7 17.6 1.71 0.186
5 Within lugar:linaje 3 19.8 6.58 0.640 0.595
6 Within Residuals 30 308. 10.3 NA NA
Valor de la tabla de F para el factor lugar con una significancia de 0.05.
qf(0.95, 1, 5)
[1] 6.607891
Valor de la tabla de F para el factor linaje con una significancia de 0.05.
qf(0.95, 3, 30)
[1] 2.922277
Conclusión.
Con respecto al Factor lugar: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor lugar tienen un efecto sobre el rendimiento estadísticamente diferente a 0.
Con respecto al Factor linaje: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor linaje tienen un efecto sobre el rendimiento estadísticamente similar a 0.
agricolae::cv.model(modelo.dbca)
[1] 23.9758
Comparaciones de medias para los efectos principales del Factor Lugar
get_df_ea(aov.dbca)
[1] 5
get_mse_ea(aov.dbca)
[1] 7.733333
data_D %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta lugar, # Cambiar según nombre de variable independienteDFerror =get_df_ea(aov.dbca), MSerror =get_mse_ea(aov.dbca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ lugar
LSD t Test for rdto
Mean Square Error: 7.733333
lugar, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
Lima 11.16667 2.443566 24 9.707486 12.62585 4 15
Pisco 15.58333 3.855168 24 14.124152 17.04251 9 26
Alpha: 0.05 ; DF Error: 5
Critical Value of t: 2.570582
least Significant Difference: 2.063594
Treatments with the same letter are not significantly different.
rdto groups
Pisco 15.58333 a
Lima 11.16667 b
Comparaciones de medias para los efectos principales del Factor Linaje
data_D %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta linaje, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dbca), MSerror = dvmisc::get_mse(modelo.dbca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ linaje
LSD t Test for rdto
Mean Square Error: 10.28333
linaje, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
1 13.66667 3.284491 12 11.776109 15.55722 9 20
2 14.66667 4.334499 12 12.776109 16.55722 10 26
3 13.41667 3.629634 12 11.526109 15.30722 9 22
4 11.75000 4.158780 12 9.859442 13.64056 4 20
Alpha: 0.05 ; DF Error: 30
Critical Value of t: 2.042272
least Significant Difference: 2.673653
Treatments with the same letter are not significantly different.
rdto groups
2 14.66667 a
1 13.66667 ab
3 13.41667 ab
4 11.75000 b
Análisis de DBCA de dos vías con análisis combinado en un sitio
Importación de datos
data_D <- data %>% dplyr::filter(lugar %in%"Lima")
Durbin-Watson test
data: modelo.dbca
DW = 2.7548, p-value = 0.1302
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
shapiro.test(rstudent(modelo.dbca))
Shapiro-Wilk normality test
data: rstudent(modelo.dbca)
W = 0.97148, p-value = 0.2889
ad.test(rstudent(modelo.dbca))
Anderson-Darling normality test
data: rstudent(modelo.dbca)
A = 0.52147, p-value = 0.1758
lillie.test(rstudent(modelo.dbca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dbca)
D = 0.10849, p-value = 0.169
Exact one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dbca)
D = 0.096609, p-value = 0.7251
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dbca))
Cramer-von Mises normality test
data: rstudent(modelo.dbca)
W = 0.089678, p-value = 0.1511
pearson.test(rstudent(modelo.dbca))
Pearson chi-square normality test
data: rstudent(modelo.dbca)
P = 9.5, p-value = 0.2187
sf.test(rstudent(modelo.dbca))
Shapiro-Francia normality test
data: rstudent(modelo.dbca)
W = 0.96843, p-value = 0.19
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dbca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.3413064, Df = 1, p = 0.55908
bptest(modelo.dbca)
studentized Breusch-Pagan test
data: modelo.dbca
BP = 24.756, df = 17, p-value = 0.1003
bptest(modelo.dbca, studentize = F)
Breusch-Pagan test
data: modelo.dbca
BP = 28.444, df = 17, p-value = 0.04001
olsrr::ols_test_breusch_pagan(modelo.dbca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
--------------------------------
Response : rdto
Variables: fitted values of rdto
Test Summary
----------------------------
DF = 1
Chi2 = 0.3413064
Prob > Chi2 = 0.5590761
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.
Estadísticas globales
modelo.dbca %>%gvlma()
Call:
lm(formula = rdto ~ epoca * linaje + bloque %in% epoca + bloque,
data = data_D)
Coefficients:
(Intercept) epoca2 linaje2 linaje3
1.108e+01 5.125e+00 -1.667e-01 -1.500e+00
linaje4 bloqueII bloqueIII bloqueIV
-3.667e+00 1.750e+00 7.500e-01 1.000e+00
bloqueV bloqueVI epoca2:linaje2 epoca2:linaje3
2.500e+00 2.500e+00 -1.167e+00 -3.833e+00
epoca2:linaje4 epoca2:bloqueII epoca2:bloqueIII epoca2:bloqueIV
-3.500e+00 1.250e+00 1.000e+00 -9.971e-16
epoca2:bloqueV epoca2:bloqueVI
7.500e-01 3.250e+00
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = .)
Value p-value Decision
Global Stat 10.1485 0.03800 Assumptions NOT satisfied!
Skewness 0.8206 0.36500 Assumptions acceptable.
Kurtosis 0.1776 0.67347 Assumptions acceptable.
Link Function 3.7952 0.05140 Assumptions acceptable.
Heteroscedasticity 5.3551 0.02066 Assumptions NOT satisfied!
\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)
\(H_1: \beta_j \neq 0\text{; en al menos un nivel del factor B.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(epoca/Bloque)” en el caso de efectos aleatorios o “Error(epoca:Bloque)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos.
Valor de la tabla de F para el factor lugar con una significancia de 0.05.
qf(0.95, 1, 5)
[1] 6.607891
Valor de la tabla de F para el factor linaje con una significancia de 0.05.
qf(0.95, 3, 30)
[1] 2.922277
Conclusión.
Con respecto al Factor epoca: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor epoca tienen un efecto sobre el rendimiento estadísticamente diferente a 0.
Con respecto al Factor linaje: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor linaje tienen un efecto sobre el rendimiento estadísticamente diferente a 0.
agricolae::cv.model(modelo.dbca)
[1] 19.33897
Comparaciones de medias para los efectos principales del Factor Lugar
get_df_ea(aov.dbca)
[1] 5
get_mse_ea(aov.dbca)
[1] 2.870833
data_D %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta epoca, # Cambiar según nombre de variable independienteDFerror =get_df_ea(aov.dbca), MSerror =get_mse_ea(aov.dbca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ epoca
LSD t Test for rdto
Mean Square Error: 2.870833
epoca, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
1 11.16667 2.443566 24 10.27761 12.05572 4 15
2 15.20833 4.242427 24 14.31928 16.09739 7 26
Alpha: 0.05 ; DF Error: 5
Critical Value of t: 2.570582
least Significant Difference: 1.257317
Treatments with the same letter are not significantly different.
rdto groups
2 15.20833 a
1 11.16667 b
Comparaciones de medias para los efectos principales del Factor Linaje
data_D %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta linaje, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dbca), MSerror = dvmisc::get_mse(modelo.dbca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ linaje
LSD t Test for rdto
Mean Square Error: 6.504167
linaje, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
1 15.58333 3.848455 12 14.079780 17.08689 11 23
2 14.83333 4.152400 12 13.329780 16.33689 10 26
3 12.16667 2.167249 12 10.663113 13.67022 9 17
4 10.16667 3.298301 12 8.663113 11.67022 4 16
Alpha: 0.05 ; DF Error: 30
Critical Value of t: 2.042272
least Significant Difference: 2.126346
Treatments with the same letter are not significantly different.
rdto groups
1 15.58333 a
2 14.83333 a
3 12.16667 b
4 10.16667 b
Análisis de DBCA de dos vías con análisis completo
Importación de datos
data_D <- data %>% dplyr::filter(lugar %in%"Lima")
Durbin-Watson test
data: modelo.dbca
DW = 2.7548, p-value = 0.1302
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
shapiro.test(rstudent(modelo.dbca))
Shapiro-Wilk normality test
data: rstudent(modelo.dbca)
W = 0.97148, p-value = 0.2889
ad.test(rstudent(modelo.dbca))
Anderson-Darling normality test
data: rstudent(modelo.dbca)
A = 0.52147, p-value = 0.1758
lillie.test(rstudent(modelo.dbca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dbca)
D = 0.10849, p-value = 0.169
Exact one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dbca)
D = 0.096609, p-value = 0.7251
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dbca))
Cramer-von Mises normality test
data: rstudent(modelo.dbca)
W = 0.089678, p-value = 0.1511
pearson.test(rstudent(modelo.dbca))
Pearson chi-square normality test
data: rstudent(modelo.dbca)
P = 9.5, p-value = 0.2187
sf.test(rstudent(modelo.dbca))
Shapiro-Francia normality test
data: rstudent(modelo.dbca)
W = 0.96843, p-value = 0.19
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dbca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.3413064, Df = 1, p = 0.55908
bptest(modelo.dbca)
studentized Breusch-Pagan test
data: modelo.dbca
BP = 24.756, df = 17, p-value = 0.1003
bptest(modelo.dbca, studentize = F)
Breusch-Pagan test
data: modelo.dbca
BP = 28.444, df = 17, p-value = 0.04001
olsrr::ols_test_breusch_pagan(modelo.dbca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
--------------------------------
Response : rdto
Variables: fitted values of rdto
Test Summary
----------------------------
DF = 1
Chi2 = 0.3413064
Prob > Chi2 = 0.5590761
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.
Estadísticas globales
modelo.dbca %>%gvlma()
Call:
lm(formula = rdto ~ epoca * linaje + bloque %in% epoca + bloque,
data = data_D)
Coefficients:
(Intercept) epoca2 linaje2 linaje3
1.108e+01 5.125e+00 -1.667e-01 -1.500e+00
linaje4 bloqueII bloqueIII bloqueIV
-3.667e+00 1.750e+00 7.500e-01 1.000e+00
bloqueV bloqueVI epoca2:linaje2 epoca2:linaje3
2.500e+00 2.500e+00 -1.167e+00 -3.833e+00
epoca2:linaje4 epoca2:bloqueII epoca2:bloqueIII epoca2:bloqueIV
-3.500e+00 1.250e+00 1.000e+00 -9.971e-16
epoca2:bloqueV epoca2:bloqueVI
7.500e-01 3.250e+00
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = .)
Value p-value Decision
Global Stat 10.1485 0.03800 Assumptions NOT satisfied!
Skewness 0.8206 0.36500 Assumptions acceptable.
Kurtosis 0.1776 0.67347 Assumptions acceptable.
Link Function 3.7952 0.05140 Assumptions acceptable.
Heteroscedasticity 5.3551 0.02066 Assumptions NOT satisfied!
\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)
\(H_1: \beta_j \neq 0\text{; en al menos un nivel del factor B.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(epoca/Bloque)” en el caso de efectos aleatorios o “Error(epoca:Bloque)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos.
Valor de la tabla de F para el factor lugar con una significancia de 0.05.
qf(0.95, 1, 5)
[1] 6.607891
Valor de la tabla de F para el factor linaje con una significancia de 0.05.
qf(0.95, 3, 30)
[1] 2.922277
Conclusión.
Con respecto al Factor epoca: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor epoca tienen un efecto sobre el rendimiento estadísticamente diferente a 0.
Con respecto al Factor linaje: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor linaje tienen un efecto sobre el rendimiento estadísticamente diferente a 0.
agricolae::cv.model(modelo.dbca)
[1] 19.33897
Comparaciones de medias para los efectos principales del Factor Lugar
get_df_ea(aov.dbca)
[1] 5
get_mse_ea(aov.dbca)
[1] 2.870833
data_D %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta epoca, # Cambiar según nombre de variable independienteDFerror =get_df_ea(aov.dbca), MSerror =get_mse_ea(aov.dbca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ epoca
LSD t Test for rdto
Mean Square Error: 2.870833
epoca, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
1 11.16667 2.443566 24 10.27761 12.05572 4 15
2 15.20833 4.242427 24 14.31928 16.09739 7 26
Alpha: 0.05 ; DF Error: 5
Critical Value of t: 2.570582
least Significant Difference: 1.257317
Treatments with the same letter are not significantly different.
rdto groups
2 15.20833 a
1 11.16667 b
Comparaciones de medias para los efectos principales del Factor Linaje
data_D %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta linaje, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dbca), MSerror = dvmisc::get_mse(modelo.dbca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ linaje
LSD t Test for rdto
Mean Square Error: 6.504167
linaje, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
1 15.58333 3.848455 12 14.079780 17.08689 11 23
2 14.83333 4.152400 12 13.329780 16.33689 10 26
3 12.16667 2.167249 12 10.663113 13.67022 9 17
4 10.16667 3.298301 12 8.663113 11.67022 4 16
Alpha: 0.05 ; DF Error: 30
Critical Value of t: 2.042272
least Significant Difference: 2.126346
Treatments with the same letter are not significantly different.
rdto groups
1 15.58333 a
2 14.83333 a
3 12.16667 b
4 10.16667 b
Análisis de DBCA de dos vías con análisis combinado en un sitio
Definición del modelo
modelo.dbca <-lm(rdto ~ epoca * linaje * lugar + bloque %in% lugar + linaje:bloque %in% epoca + bloque %in% lugar:epoca + bloque, data = data)summary(modelo.dbca)
Durbin-Watson test
data: modelo.dbca
DW = 2.3263, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
shapiro.test(rstudent(modelo.dbca))
Shapiro-Wilk normality test
data: rstudent(modelo.dbca)
W = 0.99205, p-value = 0.8432
ad.test(rstudent(modelo.dbca))
Anderson-Darling normality test
data: rstudent(modelo.dbca)
A = 0.26803, p-value = 0.6771
lillie.test(rstudent(modelo.dbca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dbca)
D = 0.065238, p-value = 0.4014
Exact one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dbca)
D = 0.058852, p-value = 0.874
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dbca))
Cramer-von Mises normality test
data: rstudent(modelo.dbca)
W = 0.04972, p-value = 0.5106
pearson.test(rstudent(modelo.dbca))
Pearson chi-square normality test
data: rstudent(modelo.dbca)
P = 13.687, p-value = 0.1877
sf.test(rstudent(modelo.dbca))
Shapiro-Francia normality test
data: rstudent(modelo.dbca)
W = 0.99385, p-value = 0.8867
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dbca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.9637412, Df = 1, p = 0.32625
bptest(modelo.dbca)
studentized Breusch-Pagan test
data: modelo.dbca
BP = 96, df = 65, p-value = 0.007467
bptest(modelo.dbca, studentize = F)
Breusch-Pagan test
data: modelo.dbca
BP = 76.736, df = 65, p-value = 0.1513
olsrr::ols_test_breusch_pagan(modelo.dbca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
--------------------------------
Response : rdto
Variables: fitted values of rdto
Test Summary
----------------------------
DF = 1
Chi2 = 0.9637412
Prob > Chi2 = 0.3262461
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.
\(Y_{ijkl}\) = Valor observado de la variable respuesta.
\(\hat{Y}_{ijkl}\) = Valor ajustado de la variable respuesta.
\(\mu\) = Promedio observado de la variable respuesta.
\(\tau_{i}\) = Efecto del i-ésimo nivel del factor epoca.
\(\beta_{j}\) = Efecto del j-ésimo nivel del factor linaje.
\(\gamma_{k}\) = Efecto del k-ésimo nivel del factor lugar.
\(\text{Error}(\gamma\text{Bloque})_{k(l)}\) = Efecto del k-ésimo nivel de lugar en el l-ésimo nivel de Bloque.
\(\text{Error}(\tau\beta\text{Bloque})_{ij(l)}\) = Efecto del i-ésimo nivel de epoca y el j-ésimo nivel de linaje en el l-ésimo nivel de Bloque.
\(\text{Error}(\tau\gamma\text{Bloque})_{ik(l)}\) = Efecto del i-ésimo nivel de epoca y el k-ésimo nivel de lugar en el l-ésimo nivel de Bloque.
\((\tau\beta)_{ij}\) = Efecto de la interacción entre el i-ésimo nivel del factor época y el j-ésimo nivel del factor linaje.
\((\tau\gamma)_{ik}\) = Efecto de la interacción entre el i-ésimo nivel del factor época y el k-ésimo nivel del factor lugar.
\((\beta\gamma)_{jk}\) = Efecto de la interacción entre el j-ésimo nivel del factor linaje y el k-ésimo nivel del factor lugar.
\((\tau\beta\gamma)_{ijk}\) = Efecto de la interacción entre el i-ésimo nivel del factor época, el j-ésimo nivel del factor linaje y el k-ésimo nivel del factor lugar
\(\delta_{l}\) = Efecto del k-ésimo nivel de Bloque.
\(\epsilon_{ijkl}\) = Residuo observado del modelo.
Pruebas de hipótesis
Para el factor A (epoca):
\(H_0: \tau_{A1} = \tau_{A2} = 0\)
\(H_1: \text{En al menos un nivel del factor A el } \tau \text{ es diferente a los demás.}\)
\(H_1: \tau_i \neq 0\text{; en al menos un nivel del factor A.}\)
\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)
\(H_1: \beta_j \neq 0\text{; en al menos un nivel del factor B.}\)
Para el factor C (lugar):
\(H_0: \gamma_{C1} = \gamma_{C2} = 0\)
\(H_1: \text{En al menos un nivel del factor C el } \gamma \text{ es diferente a los demás.}\)
\(H_1: \gamma_k \neq 0\text{; en al menos un nivel del factor C.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
Con respecto al Factor epoca: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor epoca tienen un efecto sobre el rendimiento estadísticamente diferente a 0.
Con respecto al Factor linaje: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor linaje tienen un efecto sobre el rendimiento estadísticamente diferente a 0.
Con respecto al Factor lugar: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor lugar tienen un efecto sobre el rendimiento estadísticamente diferente a 0.
agricolae::cv.model(modelo.dbca)
[1] 16.90691
Comparaciones de medias para los efectos principales del Factor Lugar
data %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta epoca, # Cambiar según nombre de variable independienteDFerror =get_df_term(aov.dbca, paste0("linaje",":","epoca")), MSerror =get_mse_term(aov.dbca, paste0("linaje",":","epoca")),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ epoca
LSD t Test for rdto
Mean Square Error: 10.49881
epoca, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
1 13.37500 3.895579 48 12.42556 14.32444 4 26
2 17.91667 4.694194 48 16.96722 18.86611 7 27
Alpha: 0.05 ; DF Error: 35
Critical Value of t: 2.030108
least Significant Difference: 1.342714
Treatments with the same letter are not significantly different.
rdto groups
2 17.91667 a
1 13.37500 b
Comparaciones de medias para los efectos principales del Factor Linaje
data %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta linaje, # Cambiar según nombre de variable independienteDFerror =get_df_term(aov.dbca, paste0("linaje",":","epoca")), MSerror =get_mse_term(aov.dbca, paste0("linaje",":","epoca")),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ linaje
LSD t Test for rdto
Mean Square Error: 10.49881
linaje, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
1 16.54167 4.736346 24 15.19895 17.88438 9 27
2 16.95833 4.685585 24 15.61562 18.30105 10 26
3 15.33333 4.410133 24 13.99062 16.67605 9 22
4 13.75000 5.219112 24 12.40729 15.09271 4 24
Alpha: 0.05 ; DF Error: 35
Critical Value of t: 2.030108
least Significant Difference: 1.898884
Treatments with the same letter are not significantly different.
rdto groups
2 16.95833 a
1 16.54167 a
3 15.33333 ab
4 13.75000 b
Comparaciones de medias para los efectos principales del Factor Lugar
data %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta lugar, # Cambiar según nombre de variable independienteDFerror =get_df_term(aov.dbca, "lugar"), MSerror =get_mse_term(aov.dbca, "lugar"),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ lugar
LSD t Test for rdto
Mean Square Error: 22.76667
lugar, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
Lima 13.18750 3.987514 48 11.41714 14.95786 4 26
Pisco 18.10417 4.415830 48 16.33381 19.87452 9 27
Alpha: 0.05 ; DF Error: 5
Critical Value of t: 2.570582
least Significant Difference: 2.503661
Treatments with the same letter are not significantly different.
rdto groups
Pisco 18.10417 a
Lima 13.18750 b
Comparaciones de medias para las interacciones linaje (B) x lugar (C)
Para los niveles del factor B dentro del nivel C1:
B1 vs B2:
\(H_0: \mu_{B1} - \mu_{B2} = 0\)
\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)
B1 vs B3:
\(H_0: \mu_{B1} - \mu_{B3} = 0\)
\(H_1: \mu_{B1} - \mu_{B3} \neq 0\)
B2 vs B3:
\(H_0: \mu_{B2} - \mu_{B3} = 0\)
\(H_1: \mu_{B2} - \mu_{B3} \neq 0\)
…
B3 vs B4:
\(H_0: \mu_{B3} - \mu_{B4} = 0\)
\(H_1: \mu_{B3} - \mu_{B4} \neq 0\)
Para los niveles del factor B dentro del nivel C1:
B1 vs B2:
\(H_0: \mu_{B1} - \mu_{B2} = 0\)
\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)
NOTA: Repetir este proceso para cada nivel de B y cada nivel de C.
multcomp.test_2factors <-function(object, respuesta, factor_name1, factor_name2, test, aov){# Función auxiliar para aplicar la prueba de comparaciones múltiples a un data frame multcomp_df <-function(df, respuesta, factor_name, test, aov){if (test =="LSD") { comp <-LSD.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[5]] } elseif (test =="HSD") { comp <-HSD.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[5]] } elseif (test =="duncan") { comp <-duncan.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[6]] } elseif (test =="SNK") { comp <-SNK.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[6]] } else {stop("Test no válido. Los tests disponibles son LSD, HSD, duncan, y SNK.") }return(comp %>%rename("y"="df[[respuesta]]") ) }# Aplicar la prueba de comparaciones múltiples a todos los data frames dentro de filters1 comp_filters1 <-lapply(object[[1]], function(df){multcomp_df(df, respuesta, factor_name2, test, aov) %>%arrange(row.names(.)) %>%rownames_to_column(var ="x") %>%relocate(x) })# Aplicar la prueba de comparaciones múltiples a todos los data frames dentro de filters2 comp_filters2 <-lapply(object[[2]], function(df){multcomp_df(df, respuesta, factor_name1, test, aov) %>% dplyr::arrange(row.names(.)) %>%rownames_to_column(var ="x") %>%relocate(x) %>% dplyr::mutate(groups =toupper(groups)) })row.names(comp_filters1) <-NULLrow.names(comp_filters2) <-NULL# Retornar una lista con las comparaciones múltiples para cada data frame result <-list() result[[as.name(substitute(factor_name1))]] <- comp_filters1 result[[as.name(substitute(factor_name2))]] <- comp_filters2return(#list(#"Comparación de los niveles del factor B dentro de cada nivel del factor A", result#[[1]],# "Comparación de los niveles del factor A dentro de cada nivel del factor B",# result[[2]]# ) )}
$linaje
$linaje$`1`
x y groups
1 Lima 15.58333 a
2 Pisco 17.50000 a
$linaje$`2`
x y groups
1 Lima 14.83333 b
2 Pisco 19.08333 a
$linaje$`3`
x y groups
1 Lima 12.16667 b
2 Pisco 18.50000 a
$linaje$`4`
x y groups
1 Lima 10.16667 b
2 Pisco 17.33333 a
$lugar
$lugar$Lima
x y groups
1 1 15.58333 A
2 2 14.83333 A
3 3 12.16667 B
4 4 10.16667 B
$lugar$Pisco
x y groups
1 1 17.50000 A
2 2 19.08333 A
3 3 18.50000 A
4 4 17.33333 A
# Convertir la lista a un data framedf <- result.comp %>% reshape2::melt() %>%ungroup() %>% dplyr::select(-variable) %>%relocate("Factor"="L1","Nivel"="L2", x,"y"="value", groups)df
Factor Nivel x y groups
1 linaje 1 Lima 15.58333 a
2 linaje 1 Pisco 17.50000 a
3 linaje 2 Lima 14.83333 b
4 linaje 2 Pisco 19.08333 a
5 linaje 3 Lima 12.16667 b
6 linaje 3 Pisco 18.50000 a
7 linaje 4 Lima 10.16667 b
8 linaje 4 Pisco 17.33333 a
9 lugar Lima 1 15.58333 A
10 lugar Lima 2 14.83333 A
11 lugar Lima 3 12.16667 B
12 lugar Lima 4 10.16667 B
13 lugar Pisco 1 17.50000 A
14 lugar Pisco 2 19.08333 A
15 lugar Pisco 3 18.50000 A
16 lugar Pisco 4 17.33333 A
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# # Convertir level1 y level2 en nombres simbólicos# level1 <- as.name(level1)# level2 <- as.name(level2)# # filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1, level2)) %>% dplyr::mutate(groups =paste0(groups.y,groups.x)) %>% dplyr::select(!!level1,!!level2, y, groups) #%>%# rename(!!level1 := x,# !!level2 := Nivel)# devolver una lista con los dos subconjuntosreturn(df)}
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1,level2)) %>% dplyr::mutate(y =paste0(round(y,2), " ", groups.y, groups.x)) %>% dplyr::select(!!level1,!!level2, y)# devolver una lista con los dos subconjuntosreturn(df)}
modelo.dca <-lm(y ~ variedad*tiempo + rep:variedad, data = data)
Nota
La expresión “rep:variedad” es para considerar la interacción como efecto aleatorio.
La expresión “rep %in% variedad” permite que se genere una estructura del análisis de varianza donde se anide el efecto de la variedad dentro de cada bloque
Durbin-Watson test
data: modelo.dca
DW = 2.0611, p-value = 0.2925
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos de la infestación de la enfermedad son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos de la infestación de la enfermedad es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos de la infestación de la enfermedad es similar a la función normal}\)
shapiro.test(rstudent(modelo.dca))
Shapiro-Wilk normality test
data: rstudent(modelo.dca)
W = 0.83592, p-value = 1.346e-09
ad.test(rstudent(modelo.dca))
Anderson-Darling normality test
data: rstudent(modelo.dca)
A = 6.6044, p-value = 2.676e-16
lillie.test(rstudent(modelo.dca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dca)
D = 0.22598, p-value = 7.152e-15
Asymptotic one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dca)
D = 0.21702, p-value = 7.634e-05
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dca))
Cramer-von Mises normality test
data: rstudent(modelo.dca)
W = 1.3725, p-value = 7.37e-10
pearson.test(rstudent(modelo.dca))
Pearson chi-square normality test
data: rstudent(modelo.dca)
P = 118.07, p-value < 2.2e-16
sf.test(rstudent(modelo.dca))
Shapiro-Francia normality test
data: rstudent(modelo.dca)
W = 0.82123, p-value = 7.281e-09
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos de la infestación de la enfermedad es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza de la infestación de la enfermedad es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza de la infestación de la enfermedad no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 130.6704, Df = 1, p = < 2.22e-16
bptest(modelo.dca)
studentized Breusch-Pagan test
data: modelo.dca
BP = 73.988, df = 47, p-value = 0.007234
bptest(modelo.dca, studentize = F)
Breusch-Pagan test
data: modelo.dca
BP = 242.17, df = 47, p-value < 2.2e-16
olsrr::ols_test_breusch_pagan(modelo.dca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
-----------------------------
Response : y
Variables: fitted values of y
Test Summary
-------------------------------
DF = 1
Chi2 = 130.6704
Prob > Chi2 = 2.923296e-30
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza de la infestación de la enfermedad es constante con respecto a los valores ajustados de la infestación de la enfermedad.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto a la infestación de la enfermedad, se debe proceder a realizar el análisis de varianza.
Estadísticas globales
modelo.dca %>%gvlma()
Call:
lm(formula = y ~ variedad * tiempo + rep:variedad, data = data)
Coefficients:
(Intercept) variedadDuraznillo
-1.7778 15.6111
variedadNegra variedadPukaÑawi
21.7778 -2.3889
variedadTumbay variedadYanaPepino
18.8889 13.4444
tiempos4 tiempos10
-9.6667 81.3333
tiempos12 tiempos15
81.0000 -12.6667
tiempos16 variedadDuraznillo:tiempos4
-15.3333 1.6667
variedadNegra:tiempos4 variedadPukaÑawi:tiempos4
2.3333 1.3333
variedadTumbay:tiempos4 variedadYanaPepino:tiempos4
6.3333 4.3333
variedadDuraznillo:tiempos10 variedadNegra:tiempos10
-96.3333 -93.6667
variedadPukaÑawi:tiempos10 variedadTumbay:tiempos10
-5.0000 -87.6667
variedadYanaPepino:tiempos10 variedadDuraznillo:tiempos12
-92.6667 -63.6667
variedadNegra:tiempos12 variedadPukaÑawi:tiempos12
-79.6667 -36.6667
variedadTumbay:tiempos12 variedadYanaPepino:tiempos12
-85.6667 -77.0000
variedadDuraznillo:tiempos15 variedadNegra:tiempos15
-5.3333 -5.0000
variedadPukaÑawi:tiempos15 variedadTumbay:tiempos15
-4.3333 -5.0000
variedadYanaPepino:tiempos15 variedadDuraznillo:tiempos16
-2.0000 -4.0000
variedadNegra:tiempos16 variedadPukaÑawi:tiempos16
-4.6667 -3.0000
variedadTumbay:tiempos16 variedadYanaPepino:tiempos16
-3.3333 0.6667
variedadCanchan:rep2 variedadDuraznillo:rep2
11.5000 3.0000
variedadNegra:rep2 variedadPukaÑawi:rep2
1.3333 47.5000
variedadTumbay:rep2 variedadYanaPepino:rep2
2.1667 0.1667
variedadCanchan:rep3 variedadDuraznillo:rep3
40.8333 13.5000
variedadNegra:rep3 variedadPukaÑawi:rep3
-1.3333 21.0000
variedadTumbay:rep3 variedadYanaPepino:rep3
5.5000 8.8333
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = .)
Value p-value Decision
Global Stat 184.8215 0.000e+00 Assumptions NOT satisfied!
Skewness 0.7003 4.027e-01 Assumptions acceptable.
Kurtosis 93.0064 0.000e+00 Assumptions NOT satisfied!
Link Function 71.1713 0.000e+00 Assumptions NOT satisfied!
Heteroscedasticity 19.9435 7.976e-06 Assumptions NOT satisfied!
\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)
\(H_1: \beta_j \neq 0\text{; en al menos un nivel del factor B.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
anova(modelo.dca, test ="F")
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
variedad 5 13725 2745.0 5.0396 0.0006392 ***
tiempo 5 27916 5583.1 10.2501 3.922e-07 ***
variedad:tiempo 25 33778 1351.1 2.4805 0.0021249 **
variedad:rep 12 13142 1095.2 2.0106 0.0387578 *
Residuals 60 32681 544.7
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Importante
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(Zona/rep)” en el caso de efectos aleatorios o “Error(Zona:rep)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos.
aov(y ~ variedad*tiempo +Error(variedad:rep), data = data) -> aov.dcasummary(aov.dca)
Error: variedad:rep
Df Sum Sq Mean Sq F value Pr(>F)
variedad 5 13725 2745 2.507 0.0892 .
Residuals 12 13142 1095
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
tiempo 5 27916 5583 10.250 3.92e-07 ***
variedad:tiempo 25 33778 1351 2.481 0.00212 **
Residuals 60 32681 545
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
broom::tidy(aov.dca)
# A tibble: 5 × 7
stratum term df sumsq meansq statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 variedad:rep variedad 5 13725. 2745. 2.51 0.0892
2 variedad:rep Residuals 12 13142. 1095. NA NA
3 Within tiempo 5 27916. 5583. 10.3 0.000000392
4 Within variedad:tiempo 25 33778. 1351. 2.48 0.00212
5 Within Residuals 60 32681. 545. NA NA
Valor de la tabla de F para el factor variedad con una significancia de 0.05.
qf(0.95, 5, 12)
[1] 3.105875
Valor de la tabla de F para el factor tiempo con una significancia de 0.05.
qf(0.95, 5, 60)
[1] 2.36827
Conclusión.
Con respecto al Factor variedad: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor Zona tienen un efecto sobre la infestación de la enfermedad estadísticamente similar a 0.
Con respecto al Factor tiempo: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor tiempo posee un efecto sobre la infestación de la enfermedad estadísticamente diferente a 0.
Con respecto a la interacción variedad*tiempo: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos una interacción entre variedad y tiempo tuvo un efecto de antagonismo o sinergismo en infestación de la enfermedad.
agricolae::cv.model(modelo.dca)
[1] 127.6884
Comparaciones de medias para los efectos principales del Factor variedad
data %>%with(LSD.test( y, # Cambiar según nombre de variable respuesta variedad, # Cambiar según nombre de variable independienteDFerror =get_df_ea(aov.dca), MSerror =get_mse_ea(aov.dca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: y ~ variedad
LSD t Test for y
Mean Square Error: 1095.157
variedad, means and individual ( 95 %) CI
y std r LCL UCL Min Max
Canchan 36.444444 53.575919 18 19.449414 53.43947 0 178
Duraznillo 12.166667 19.147262 18 -4.828364 29.16170 0 81
Negra 10.666667 8.977095 18 -6.328364 27.66170 0 30
PukaÑawi 31.500000 53.140491 18 14.504970 48.49503 0 157
Tumbay 11.222222 8.142521 18 -5.772808 28.21725 0 23
YanaPepino 7.666667 10.803050 18 -9.328364 24.66170 0 43
Alpha: 0.05 ; DF Error: 12
Critical Value of t: 2.178813
least Significant Difference: 24.0346
Treatments with the same letter are not significantly different.
y groups
Canchan 36.444444 a
PukaÑawi 31.500000 ab
Duraznillo 12.166667 b
Tumbay 11.222222 b
Negra 10.666667 b
YanaPepino 7.666667 b
Comparaciones de medias para los efectos principales del Factor tiempo
data %>%with(LSD.test( y, # Cambiar según nombre de variable respuesta tiempo, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dca), MSerror = dvmisc::get_mse(modelo.dca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: y ~ tiempo
LSD t Test for y
Mean Square Error: 544.6907
tiempo, means and individual ( 95 %) CI
y std r LCL UCL Min Max
s10 36.7777778 56.8667747 18 25.774212964 47.78134 0 178
s12 41.8888889 47.1142402 18 30.885324075 52.89245 1 157
s15 1.7222222 2.7824144 18 -9.281342591 12.72579 0 9
s16 0.2777778 0.7519039 18 -10.725787036 11.28134 0 3
s3 18.0000000 2.8697202 18 6.996435187 29.00356 10 20
s4 11.0000000 4.4325003 18 -0.003564813 22.00356 3 18
Alpha: 0.05 ; DF Error: 60
Critical Value of t: 2.000298
least Significant Difference: 15.56139
Treatments with the same letter are not significantly different.
y groups
s12 41.8888889 a
s10 36.7777778 a
s3 18.0000000 b
s4 11.0000000 bc
s15 1.7222222 c
s16 0.2777778 c
Comparaciones de medias para las interacciones
Para los niveles del factor A dentro del nivel B1:
A1 vs A2:
\(H_0: \mu_{A1} - \mu_{A2} = 0\)
\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)
…
A5 vs A6:
\(H_0: \mu_{A5} - \mu_{A6} = 0\)
\(H_1: \mu_{A5} - \mu_{A6} \neq 0\)
Para los niveles del factor B dentro del nivel A1:
B1 vs B2:
\(H_0: \mu_{B1} - \mu_{B2} = 0\)
\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)
…
B5 vs B6:
\(H_0: \mu_{B5} - \mu_{B6} = 0\)
\(H_1: \mu_{B5} - \mu_{B6} \neq 0\)
NOTA: Repetir este proceso para cada nivel de A y cada nivel de B.
Análisis de varianza para interacción de dos factores con el paquete phia
Comparación de los niveles de B dentro de cada nivel de A
multcomp.test_2factors <-function(object, respuesta, factor_name1, factor_name2, test, aov, reverse){# Función auxiliar para aplicar la prueba de comparaciones múltiples a un data frame multcomp_df <-function(df, respuesta, factor_name, test, aov){if (test =="LSD") { comp <-LSD.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[5]] } elseif (test =="HSD") { comp <-HSD.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[5]] } elseif (test =="duncan") { comp <-duncan.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[6]] } elseif (test =="SNK") { comp <-SNK.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[6]] } else {stop("Test no válido. Los tests disponibles son LSD, HSD, duncan, y SNK.") }if (reverse ==TRUE) { comp <- comp %>% dplyr::mutate(groups =rev(groups)) } elseif (reverse ==FALSE) { comp <- comp } else {stop("Test no válido. Los tests disponibles son LSD, HSD, duncan, y SNK.") }return(comp %>%rename("y"="df[[respuesta]]") ) }# Aplicar la prueba de comparaciones múltiples a todos los data frames dentro de filters1 comp_filters1 <-lapply(object[[1]], function(df){multcomp_df(df, respuesta, factor_name2, test, aov) %>%arrange(row.names(.)) %>%rownames_to_column(var ="x") %>%relocate(x) })# Aplicar la prueba de comparaciones múltiples a todos los data frames dentro de filters2 comp_filters2 <-lapply(object[[2]], function(df){multcomp_df(df, respuesta, factor_name1, test, aov) %>% dplyr::arrange(row.names(.)) %>%rownames_to_column(var ="x") %>%relocate(x) %>% dplyr::mutate(groups =toupper(groups)) })row.names(comp_filters1) <-NULLrow.names(comp_filters2) <-NULL# Retornar una lista con las comparaciones múltiples para cada data frame result <-list() result[[as.name(substitute(factor_name1))]] <- comp_filters1 result[[as.name(substitute(factor_name2))]] <- comp_filters2return(#list(#"Comparación de los niveles del factor B dentro de cada nivel del factor A", result#[[1]],# "Comparación de los niveles del factor A dentro de cada nivel del factor B",# result[[2]]# ) )}
$variedad
$variedad$Canchan
x y groups
1 s10 97.0000000 b
2 s12 96.6666667 b
3 s15 3.0000000 a
4 s16 0.3333333 a
5 s3 15.6666667 b
6 s4 6.0000000 b
$variedad$Duraznillo
x y groups
1 s10 4.333333 a
2 s12 36.666667 a
3 s15 1.333333 a
4 s16 0.000000 a
5 s3 19.333333 a
6 s4 11.333333 a
$variedad$Negra
x y groups
1 s10 7.666667 a
2 s12 21.333333 a
3 s15 2.333333 a
4 s16 0.000000 a
5 s3 20.000000 a
6 s4 12.666667 a
$variedad$PukaÑawi
x y groups
1 s10 95.0000000 b
2 s12 63.0000000 b
3 s15 1.6666667 a
4 s16 0.3333333 a
5 s3 18.6666667 b
6 s4 10.3333333 b
$variedad$Tumbay
x y groups
1 s10 13.33333 a
2 s12 15.00000 a
3 s15 2.00000 a
4 s16 1.00000 a
5 s3 19.66667 a
6 s4 16.33333 a
$variedad$YanaPepino
x y groups
1 s10 3.333333 a
2 s12 18.666667 a
3 s15 0.000000 a
4 s16 0.000000 a
5 s3 14.666667 a
6 s4 9.333333 a
$tiempo
$tiempo$s3
x y groups
1 Canchan 15.66667 A
2 Duraznillo 19.33333 A
3 Negra 20.00000 A
4 PukaÑawi 18.66667 A
5 Tumbay 19.66667 A
6 YanaPepino 14.66667 A
$tiempo$s4
x y groups
1 Canchan 6.000000 A
2 Duraznillo 11.333333 A
3 Negra 12.666667 A
4 PukaÑawi 10.333333 A
5 Tumbay 16.333333 A
6 YanaPepino 9.333333 A
$tiempo$s10
x y groups
1 Canchan 97.000000 B
2 Duraznillo 4.333333 A
3 Negra 7.666667 B
4 PukaÑawi 95.000000 B
5 Tumbay 13.333333 B
6 YanaPepino 3.333333 A
$tiempo$s12
x y groups
1 Canchan 96.66667 C
2 Duraznillo 36.66667 C
3 Negra 21.33333 BC
4 PukaÑawi 63.00000 C
5 Tumbay 15.00000 A
6 YanaPepino 18.66667 AB
$tiempo$s15
x y groups
1 Canchan 3.000000 A
2 Duraznillo 1.333333 A
3 Negra 2.333333 A
4 PukaÑawi 1.666667 A
5 Tumbay 2.000000 A
6 YanaPepino 0.000000 A
$tiempo$s16
x y groups
1 Canchan 0.3333333 A
2 Duraznillo 0.0000000 A
3 Negra 0.0000000 A
4 PukaÑawi 0.3333333 A
5 Tumbay 1.0000000 A
6 YanaPepino 0.0000000 A
# Convertir la lista a un data framedf <- result.comp %>% reshape2::melt() %>%ungroup() %>% dplyr::select(-variable) %>%relocate("Factor"="L1","Nivel"="L2", x,"y"="value", groups) #%>%# group_by(Factor, Nivel) %>%# dplyr::mutate(groups = rev(groups)) %>%# dplyr::ungroup()df
Factor Nivel x y groups
1 variedad Canchan s10 97.0000000 b
2 variedad Canchan s12 96.6666667 b
3 variedad Canchan s15 3.0000000 a
4 variedad Canchan s16 0.3333333 a
5 variedad Canchan s3 15.6666667 b
6 variedad Canchan s4 6.0000000 b
7 variedad Duraznillo s10 4.3333333 a
8 variedad Duraznillo s12 36.6666667 a
9 variedad Duraznillo s15 1.3333333 a
10 variedad Duraznillo s16 0.0000000 a
11 variedad Duraznillo s3 19.3333333 a
12 variedad Duraznillo s4 11.3333333 a
13 variedad Negra s10 7.6666667 a
14 variedad Negra s12 21.3333333 a
15 variedad Negra s15 2.3333333 a
16 variedad Negra s16 0.0000000 a
17 variedad Negra s3 20.0000000 a
18 variedad Negra s4 12.6666667 a
19 variedad PukaÑawi s10 95.0000000 b
20 variedad PukaÑawi s12 63.0000000 b
21 variedad PukaÑawi s15 1.6666667 a
22 variedad PukaÑawi s16 0.3333333 a
23 variedad PukaÑawi s3 18.6666667 b
24 variedad PukaÑawi s4 10.3333333 b
25 variedad Tumbay s10 13.3333333 a
26 variedad Tumbay s12 15.0000000 a
27 variedad Tumbay s15 2.0000000 a
28 variedad Tumbay s16 1.0000000 a
29 variedad Tumbay s3 19.6666667 a
30 variedad Tumbay s4 16.3333333 a
31 variedad YanaPepino s10 3.3333333 a
32 variedad YanaPepino s12 18.6666667 a
33 variedad YanaPepino s15 0.0000000 a
34 variedad YanaPepino s16 0.0000000 a
35 variedad YanaPepino s3 14.6666667 a
36 variedad YanaPepino s4 9.3333333 a
37 tiempo s3 Canchan 15.6666667 A
38 tiempo s3 Duraznillo 19.3333333 A
39 tiempo s3 Negra 20.0000000 A
40 tiempo s3 PukaÑawi 18.6666667 A
41 tiempo s3 Tumbay 19.6666667 A
42 tiempo s3 YanaPepino 14.6666667 A
43 tiempo s4 Canchan 6.0000000 A
44 tiempo s4 Duraznillo 11.3333333 A
45 tiempo s4 Negra 12.6666667 A
46 tiempo s4 PukaÑawi 10.3333333 A
47 tiempo s4 Tumbay 16.3333333 A
48 tiempo s4 YanaPepino 9.3333333 A
49 tiempo s10 Canchan 97.0000000 B
50 tiempo s10 Duraznillo 4.3333333 A
51 tiempo s10 Negra 7.6666667 B
52 tiempo s10 PukaÑawi 95.0000000 B
53 tiempo s10 Tumbay 13.3333333 B
54 tiempo s10 YanaPepino 3.3333333 A
55 tiempo s12 Canchan 96.6666667 C
56 tiempo s12 Duraznillo 36.6666667 C
57 tiempo s12 Negra 21.3333333 BC
58 tiempo s12 PukaÑawi 63.0000000 C
59 tiempo s12 Tumbay 15.0000000 A
60 tiempo s12 YanaPepino 18.6666667 AB
61 tiempo s15 Canchan 3.0000000 A
62 tiempo s15 Duraznillo 1.3333333 A
63 tiempo s15 Negra 2.3333333 A
64 tiempo s15 PukaÑawi 1.6666667 A
65 tiempo s15 Tumbay 2.0000000 A
66 tiempo s15 YanaPepino 0.0000000 A
67 tiempo s16 Canchan 0.3333333 A
68 tiempo s16 Duraznillo 0.0000000 A
69 tiempo s16 Negra 0.0000000 A
70 tiempo s16 PukaÑawi 0.3333333 A
71 tiempo s16 Tumbay 1.0000000 A
72 tiempo s16 YanaPepino 0.0000000 A
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# # Convertir level1 y level2 en nombres simbólicos# level1 <- as.name(level1)# level2 <- as.name(level2)# # filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1, level2)) %>% dplyr::mutate(groups =paste0(groups.y,groups.x)) %>% dplyr::select(!!level1,!!level2, y, groups) #%>%# rename(!!level1 := x,# !!level2 := Nivel)# devolver una lista con los dos subconjuntosreturn(df)}
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1,level2)) %>% dplyr::mutate(y =paste0(round(y,2), " ", groups.y, groups.x)) %>% dplyr::select(!!level1,!!level2, y)# devolver una lista con los dos subconjuntosreturn(df)}
